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Preface 


The book is intended to provide an introduction to the effects of nonlinear elements in feedback control systems. A 
central topic is the use of the Describing Function (DF) method since in combination with simulation it provides an 
excellent approach for the practicing engineer and follows on logically from a first course in classical control, such as 
the companion volume in this series. Some of the basic material on the topic can be found in my earlier book which is 


frequently referenced throughout the text. 


The first chapter provides an introduction to nonlinearity from the basic definition to a discussion of the possible effects 
it can have on a system and the different behaviour that might be found, in particular when it occurs within a cascade 


system of elements or in a feedback loop. The final part of the chapter gives a brief overview of the contents of the book. 


Phase plane methods for second order systems are covered in the second chapter. Many systems, particularly 
electromechanical ones, can be approximated by second order models so the concept can be particularly useful in practice. 
The second order linear system is first considered as surprisingly it is rarely covered in linear control system texts. The 
method has the big advantage in that the effects of more than one nonlinear element may be considered. The study is 


supported by several simulations done in Simulink including one where sliding motion takes place. 


The third chapter is the first of three devoted to the study of feedback loops using DF methods. Although the method is 
an approximate technique its value, and limitations, are supported by a large number of examples containing analytical 
results and simulations, including the estimation of limit cycles and loop stability. Many early papers on DFs showing how 
theories could be used to predict specific phenomena were supported by simulations done on analogue computers. Here 
results from digital simulations using Simulink are presented and this allows much more control of initial conditions to 
show how different modes may exist dependent on the initial conditions. In particular, Chapter 5 contains some more 
advanced work including some new results on jump resonance, so some readers may wish to omit this chapter on a first 


reading. 


A relay is a unique nonlinear element in that its output does not depend upon the input at all times but is determined by 
when the input passes through the relay switching levels. It is this feature which allows the exact determination of limit 
cycles and their stability in a feedback loop. The basic theory is presented and some simple examples covered. It is then 
shown in section 6.6 how the approach can be used with computational support to analyse quite complicated periodic 
modes in relay systems. Among the topics covered for the first time in a textbook are the evaluation of limit cycles with 
multiple pulses per cycle, as found in a satellite attitude control system, the determination of a limit cycle with sliding 


and other more advanced aspects which some readers may wish to omit. 


A practical method developed in recent years for finding suitable parameters, or tuning, for a controller based on the 
so called loop cycling method of Ziegler and Nichols has used relay produced limit cycles. Chapter 7 covers these ideas 


using both approximate analysis based on the DF and exact analysis based on the relay methods of the previous chapter. 
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Chapter 8 covers the topic of absolute stability namely trying to obtain necessary and sufficient conditions for the stability 
of a feedback loop with a single nonlinear element. This problems has exercised the minds of theorists for nearly a century 
but a solution seems no nearer! Several necessary but not sufficient results are presented, which dependent on ones 
viewpoint may be regarded as ‘conservative’ or ‘robust. The former is usually the case when one has a mathematically 
defined nonlinearity and the latter may be used because the result gives stability for any nonlinearity with certain properties, 


for example, lying within a sector. 


The final chapter, chapter 9, discusses quite briefly various methods which can be used for the design of nonlinear systems. 
The intent has been to provide sufficient information on the methods and their possible advantages and disadvantages. 
Several of them have complete books written on the topic and more detail could not have been given without the coverage 


of more specialised mathematics. 


Finally my thanks to the University of Sussex for the use of an office and access to the computing facilities during my 
retirement; to my good friend Keith Godfrey and his student Roland for the computational input on jump resonance and 


their companionship at the races and to my wife, Constance, for her love and support. 


Derek P Atherton 
University of Sussex 
Brighton 

May 2011. 
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1 Introduction 


1.1. What is nonlinearity? 


In order to analyse the behaviour of engineering systems mathematical models are required for the various components. It 
is common practice to try and obtain linear models as a rich mathematical theory exists for linear systems. These models 
will always be approximate, although possibly quite accurate for defined ranges of system variables, but inevitably nonlinear 
effects will eventually be found for large excursions of system variables. Linear systems have the important property 
that they satisfy the superposition principle. This leads to many important advantages in methods for their analysis. For 
example, in circuit theory when an RLC circuit has both d.c and a.c input voltages, the voltage or current elsewhere in 
the circuit can be found by summing the results of separate analyses for the d.c. and a.c. inputs taken individually, also if 


the magnitude of the a.c. voltage is doubled then the a.c. voltages and currents elsewhere in the circuit will be doubled. 


Thus, mathematically a linear system with input x(f) and output y(f) satisfies the property that the output for an 
input ax: (t) + bx.(t) is ay: (t) + by2(t), if y:(f) and y2(f) are the outputs in response to the inputs x: (f) and x(t) 
, respectively, and a and b are constants. A further point is that when the a.c. voltage is sinusoidal, which is normally 
assumed in using the term a.c., then all the other voltages and currents in the circuit are sinusoidal with the same frequency 


as the input. 


A nonlinear system is defined as one which does not satisfy the superposition property. The simplest form of nonlinear 
system is the static nonlinearity where the output depends only on the current value of input but in a nonlinear manner, 


for example 
y(t) = cx(t) + ax (t) (1.1) 


where the output is the summation of a linear and a nonlinear, cubic, term. This description of nonlinearity is given by a 
simple mathematical expression. Practical nonlinearities which occur in control engineering, however, can often not be 
so easily described; so that approximating them may require more complex mathematical models, such as a high order 


power series or a linear segmented approximation. 
More commonly a nonlinear differential equation, for example 
d’y(t)/d(t)’ + c(dy(t)/dt)* + y(t) = x(t) (1.2) 


will describe the behavior. From an engineering viewpoint it may be desirable to think of this equation in terms of a 
block diagram, as shown in Figure 1.1, consisting of linear dynamic elements and a static nonlinearity, in this case a 
cubic, which with input dy(t)/dt gives an output c (dy(t)/ dt)’. Unfortunately there is no general approach to solving 
nonlinear, unlike linear, differential equations. The major point about nonlinear systems, however, is that their response 
is amplitude dependent so that if a particular form of response, or some measure of it, occurs for one input magnitude 


it may not result for some other input magnitude. 
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Figure 1.1 Block diagram in Simulink for equation 1.2 


A further and very important point, is that unlike a linear operation, a nonlinear operation on a sinusoid of frequency, 
f will not produce an output at frequency f, alone. For example, if such a sinusoid is applied to the static nonlinearity of 
equation (1.1) itis easy to show from substituting x(t) = a cos(wt) where @ = 27 that there are outputs at frequencies 
fand 3f of magnitudes ca + 3da°/a and da’ /4, respectively. Perhaps the most interesting aspect of nonlinear systems 


is that, as will be shown later, they exhibit several forms of unique behaviour which are not possible in linear systems. 


1.2 Forms of nonlinearity 


All practical systems are nonlinear and in this section a brief overview is given of some nonlinear effects that often occur. 
In initial designs it may be possible to approximate the nonlinear effects by linear models but invariably it will be necessary 
to finally check their effects on the system performance either in simulation or/and the real hardware. Today’s digital 
simulation languages are very good but to use them efficiently for investigating the effects of nonlinearity, or to assist in 


the design of a nonlinear system, requires a knowledge of the supporting theoretical methods presented in this book. 


In typical control engineering problems nonlinearity may occur in the dynamics of the plant to be controlled or in the 
components used to implement the control. In the latter case, for example, a valve actuator may have a dead zone due to 
friction effects and will certainly saturate for large inputs, so this may be referred to as an inherent nonlinearity, because 
it exists although one might possibly prefer this not to be the case. Alternatively one may have intentional nonlinearities 
which have been purposely designed into the system to improve the system specifications, either for technical or economic 
reasons. A good example of this is the on-off control used in many temperature control systems, where the objective is 


to have the temperature oscillate about the required value. 


Identifying the precise form of a nonlinearity may not be easy and like all modeling exercises the golden rule is to be aware 
of the approximations in a nonlinear model and the conditions for its validity. It might be argued that linear systems theory 
is not applicable to practical control engineering problems because they are always nonlinear. This is an overstatement, 
of course, but a valid reminder. All systems have actuator saturation and in some cases it might occur for relatively low 
error signals, for example in rotary position control it is not unusual for a step input of say 10°, or even less, to produce 
the maximum motor drive torque. It simply is the result of good economical design as to produce linear operation to a 
higher torque level would require a larger and more expensive motor. Valves used to control fluid or gas flow, apart from 
having the nonlinear effects mentioned above, can have a slightly different behavior when opening compared with closing 


due to the unidirectional pressure of the fluid. 
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Friction always occurs in mechanical systems and is very difficult to model, with many quite sophisticated models having 
been presented in the literature. The simplest is to assume the three components, illustrated in Figure 1.2, of stiction, an 
abbreviation for static friction, Coulomb friction and viscous friction. As its name implies stiction is assumed to exist only 
at zero differential speed between the two contact surfaces. Coulomb friction with a value less than stiction is assumed 
to be constant at all speeds, and viscous friction is a linear effect being directly proportional to speed. In practice there 
is often a term proportional to a higher power of speed, and this is also the situation for many shaft loads, for example a 


fan for which the drive torque typically increases as a power of speed. 


Friction Force 


Viscous 


Stiction 


Coulomb 


Speed 


Figure 1.2 Three basic friction components 


A mathematical expression sometimes used to approximate friction is 
n(v) = (a— be “+ d\ vl\)sgn(v) (3) 


Here the parameters a, b, c and d are chosen to provide the desired shape against speed which initially decreases from 


the value at zero speed and then starts to increase. 


Many mechanical loads are driven through gearing rather than directly. Although geared drives, like all areas of technology, 
have improved through the years they always have some small backlash. This may be avoided by using anti-backlash 
gears, which are only available for low torques. Backlash, which is a very complicated phenomenon, involving impacts 
between surfaces is often modeled in a very simplistic manner. For example, the simple approach used in some digital 
simulation languages, such as Simulink, the simulation component of the well-known software package, Matlab, consists 


of an input-output position characteristic of two parallel straight lines with possible horizontal movement between them. 


This makes two major assumptions, first that the load shaft friction is high enough for contact to be maintained with the 
drive side of the backlash when the drive slows down to rest. Secondly when the drive reverses the backlash is crossed 
and the new drive side of the gear ‘picks up’ the load instantaneously with no loss of energy in the impact and both then 
move at the drive shaft speed. Clearly both these assumptions are never true in practice but no checks exist in Simulink 
to determine how good they are or, indeed, when they are completely invalid. Better facilities for modeling phenomena 
such as backlash can be found in simulation languages such as 20-Sim or Dynast, which do not use the block concept of 


Simulink. A good discussion of friction and backlash can be found in referencel.1. 
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Probably, the most widely used intentional nonlinearity is the relay. The on-off type, which can be described mathematically 
by the signum function, that is switches on if its input exceeds zero and off if it goes below zero, is widely used normally 
with some hysteresis between the switching levels. Use of this approach provides a control strategy where the controlled 
variable oscillates about the desired level. The switching mechanism varies significantly according to the application from 
electromechanical relays at low speed to fast electronic switches employing transistors or thyristors. A common usage of 
the relay is in the temperature control of buildings, where typically the switching is provided from a temperature sensor 
having a pool of mercury on a metal expansion coil. As the temperature drops the coil contracts and this causes a change 
in angle of the mercury capsule so that eventually the mercury moves and closes a contact. When the temperature increases 


the coil expands causing a change in angle for the mercury to flow and break the contact. 


Electronic switching controllers are being used in many modern electric motor drive systems, for example, to regulate 
phase currents in stepping motors and switched reluctance motors and to control currents in vector control drives for 
induction motors. Relays with a dead zone, that is, three position relays giving positive, negative and a zero output are 
also used. When used in a position control system the zero output allows for a steady state position within the dead zone 


but this affects the resulting steady state control accuracy. 


Hysteresis effects in magnetic materials sometimes have to be modeled. This is often done by assuming a hysteresis loop 
of the form of a B-H loop for a magnetic material typically obtained for a sinusoidal input. However the shape usually 
varies with the amplitude and frequency of the input and does not in fact remain constant with a random excitation. 
Provided the input is sinusoidal and the shape does remain reasonably constant then a nonlinear function of the form 


n(x,x) may provide a reasonable model as the path taken around the hysteresis depends upon whether the input, x, is 


increasing or decreasing, i.e the derivative of x, x. 
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1.3. Structure and Behaviour 


From a control engineering viewpoint there are two major reasons why one needs to know about nonlinearity. Firstly 
with respect to obtaining mathematical models of devices, particularly if identification techniques are being used, and 
secondly for ensuring the design meets the desired specifications when the control system is nonlinear. To appreciate 
these aspects it is appropriate to discuss very briefly in the introduction a few aspects of behaviour due to the presence 
of nonlinearity. These are dependent on the structure of the nonlinear system and the relevance of this is explained by 


again considering a sinusoidal input. 


Consider a nonlinearity x + dx’, and a dynamic element with transfer function K/s(s + b), where b and d are constants. 
If they are placed in cascade and a sinusoidal input applied the output will be a deterministic waveform containing 
two frequency components one at the same frequency as the input and the other at three times that frequency. Any 
cascade combination of linear and nonlinear elements will always produce a deterministic output for any given discrete 
input frequency spectrum, which in principle can be evaluated. New frequency components can only be created by the 


nonlinearities and the linear elements simple alter the relative magnitudes and phases of these components. 


For example, if the above nonlinearity is placed before and after the linear transfer function and a sinusoid of frequency, 
fis applied at the input, then the input to the second nonlinearity will consist of the fundamental, f, plus third harmonic, 
3f, with magnitudes and phases dependent on both the input sinusoidal magnitude and frequency. These two frequencies 
applied to the second nonlinearity will produce an output containing the frequency components, f, 3f 5f 7f, and 9f. One 
could define a frequency response for such a cascade structure of linear and nonlinear elements as the ratio of the output 
at the fundamental frequency, f, to the input sinusoid at this frequency. The result, as for a linear system, would be a 
magnitude and phase, which varies with, f, but because of the nonlinearities it would also vary with the amplitude of the 
input sinusoid. Thus an approximate frequency response model for the combination could be portrayed graphically by 
a set of frequency response plots for different input amplitudes, or gain and phase plots against amplitude for different 


frequencies. 


For many problems encountered in control engineering this may prove to be a reasonably good model since many of 
the linear dynamic elements, like the one given, have low pass dynamics so that the frequency, f, will predominate at the 
output. With no linear dynamic elements in the combination then these latter plots would be the same for all frequencies 
and the approximate, first harmonic, or quasi-linearized, model would be gain and phase curves as a function of the input 
amplitude. This representation of a nonlinear element is known as a describing function (DF), which is the topic of several 
chapters. Sometimes the above model for the combination based on the fundamental frequency, f, only is referred to as 


an amplitude and frequency dependent describing function. 
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If alternatively we assume the system structure to consist of a feedback combination of the linear and nonlinear elements, 
with the transfer function K/s(s + b) in the forward loop and the feedback loop containing x + dx° fed back through a 
negative gain then, a very different situation is possible for the response to a sinusoidal input. Dependent on the values 
of b, d and the amplitude and frequency of the input, some possibilities for the output are that it is (a) approximately 
sinusoidal with the same frequency as the input, similar to the aforementioned cascade connection; (b) approximately 
sinusoidal with a frequency related to that of the input; (c) a combination of primarily the input frequency and another 
frequency or (d) a waveform known as chaotic, which is not definable mathematically but completely repeatable for the 


same initial conditions. 


These behaviours and others are unique to nonlinear feedback systems, aspects which make such systems extremely 
interesting. However, it has meant that no general analytical method is available for predicting their behavior. Several 
approaches will be considered in this book all of which will be restricted in their applicability or, put alternatively, the 
situations which they can address. Thus the importance of simulation studies for investigating nonlinear systems in 
association with analytical methods cannot be underestimated. Much of the support for the theoretical material presented 
in the early chapters, particularly in the 50s to 60s, was done using analogue simulation. Today simulations are done 
digitally and several are included using Simulink in the following chapters to illustrate the concepts and provide solutions 
for specific problems. Care has to be taken in simulating nonlinear systems particularly those with linear segmented 


characteristics because of the discontinuities. Some comments are made on the simulations where appropriate. 


1.4 Overview of contents 


The phase plane approach discussed in chapter 2 is very useful for step response and stability studies but is basically 
restricted to second order systems. However, many engineering systems, particularly in the mechatronics field, may be 
approximated by a second order differential equation so the results are still of value. It also provides a simple basis for 


understanding some of the more advanced topics, such as optimum control and sliding mode control, covered in chapter 9. 


Stability of a feedback loop is of major importance so that it is not surprising that much early work was concerned with 
this topic. During the 1940’s engineers in several countries developed what has become known as the describing function 
method where a nonlinearity is replaced by an amplitude dependent gain to a sinusoid, known as the describing function, 
DE Chapter 3 introduces the DF, shows how its value can be calculated for various nonlinearities and includes a table of 
results. The next two chapters, 4 and 5, deal with applications of the describing function for estimating limit cycles and 
the stability of a nonlinear feedback loop. It is also shown how the describing function can be evaluated for other than 
a single sinusoidal input and applications of describing functions for a bias plus sinusoid and two sinusoids are given. 


These include areas such as limit cycle stability, jump resonance and subharmonic oscillations. 


The relay is a special form of nonlinear element where the output is not continuously dependent on the input. This allows 
special results to be developed for relay systems a topic covered in chapter 6. Chapter 7 deals with a design approach that 
has received much attention in recent years, namely using the information contained in relay produced limit cycles to set 
the parameters of a controller, typically known as relay autotuning. The topic of absolute stability, namely the development 
of exact criteria for guaranteeing the stability of a feedback loop with a single nonlinearity and linear transfer function, 
is covered in chapter 8. Dependent on the viewpoint these results may be regarded as robust, as they prove stability for 


variations in the nonlinearity, or conservative if one is considering stability for a specific nonlinearity. 
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The coverage to this point, apart from chapter 7, has like many textbooks on linear control, been primarily concerned 
with introducing analytical tools. Chapter 9, however, focuses more on design and looks at how some of the ideas covered 
and other techniques may be used for nonlinear control system design. The coverage of these topics is quite brief, some 
necessarily so because of the additional theoretical concepts which would have to be introduced to go into them more 
deeply. Hopefully, sufficient information, together with the references, is given for the reader to understand the concepts 


involved, their possible relevance for particular applications and how they might be applied. 


1.5 References 


1.1 Friedland, B, 1996 Advanced control system design, Chapter 7 Prentice Hall. 
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2 The Phase Plane Method 


2.1. Introduction 


A significant amount of the research into nonlinear differential equations in the nineteenth and early twentieth centuries 
done by mathematicians and physicists was devoted to second order differential equations. There were two major reasons 
for this, namely that the dynamics of many problems of practical interest could be approximated by these equations and 
secondly the phase plane approach allowed a graphical examination of their solutions. The systems of interest in the late 
nineteenth and early twentieth centuries were found in fields such as celestial mechanics, nonlinear mechanical systems 
and electronic oscillations. This section will introduce the basic concepts of the phase plane approach and then give a 


brief overview of how the method has been further developed for use in control system analysis and design. 


A significant amount of the early development in control theory from 1930 was driven by two areas, namely to achieve 
better control of industrial processes and to achieve better performance in fire control problems. The problems in the latter 
area received significantly more attention in the war years after 1939, where the requirement in many cases was related 
to the position control of radar antennas and guns in both stationary and moving situations. The dynamic equations 
representing many of these position control systems could be represented reasonably accurately by second order nonlinear 
differential equations. It was therefore not surprising that much of the early work on nonlinear control used the phase 
plane approach. Control engineers did make significant contributions to this field since, whereas the earlier work had 
typically assumed nonlinearities defined by continuous mathematical functions, for control system analysis it was often 
more appropriate to approximate intrinsic nonlinearity, such as friction, or intentionally introduced nonlinearity, such as a 
relay, by linear segmented characteristics. The approach is still useful today because of the fundamental insight it provides 
into aspects of nonlinear system behaviour; the fact there are still many control problems which can be approximated by 


second order dynamics, and also because more than one nonlinearity can be considered. 


2.2 Basic Principles 


The formulation used in early work on second order systems was to assume a representation in terms of the two first 


order equations 


x = P(x, x2) (2.1) 
x = O(x, x2) 
Equilibrium, or singular, points represent a stationary system for the dynamics and occur when x; = x2 = 0 
The slope of any solution curve, or trajectory, in the x, - x, state plane is 
' X14 X: 
dx. _ x2 _ Q(xi,2) (2.2) 


dx, x1 ~ P(x,X2) 


Before considering the nonlinear case further it is important to first examine the linear case, which is rarely discussed in 
detail in texts on linear control. Most, usually, just restrict the contribution to considering the step response of a second 


order system. 
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2.2.1 The Linear Case 


The equations for the linear situation may be written 


x1 = ax, + bx2 


(2.3) 
x2 = cx, + dxo 
or in matrix form 
x = Ax 
ab (2.4) 
where A = : 
cd 


Eliminating x2 from equations (2.3) by differentiating the first equation and substituting for x. from the second equation 


and x2 from the first yields 

x, — (a+ d)x + (ad — bc)x; = 0 

which after taking Laplace transforms has the characteristic equation 
s —(at+d)s+ad—be = 0) 


which is the same as the eigenvalue equation of the matrix A, with ) typically replacing s when discussed in mathematical 


texts. 

The roots, or eigenvalues, of the equation depend on the values of the elements of A. If the equation is written in the form 
V+ BA+C=0 

then the two eigenvalues are given by 

A, =— B+ (BY — 40)'" 12. and a, =— B - (B’ — 40)" /2 


The solution to the differential equation is of the form 
XI (ft) _ Ke" + Kye?! for A # Ad 


XI (t) = ke + tKye* for A — Ar 


For the linear system there is only one singular point at the origin x, = 0, x2 = 0 of the state plane and the behaviour 


of the motion near to the singular point depends on the eigenvalues. 
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Trajectories cannot intersect and can only meet at a singular point. The following four cases can occur for the singular point. 


(i) A; and A: are both real and have the same sign (either positive or negative). 
(ii) A; and A are both real and have opposite signs. 
(iii) Ai and Az are complex conjugates with non zero real parts. 


(iv) A; and A> are imaginary, that is complex conjugates with zero real parts. 


Case (i) corresponds to a singular point known as a node, which is stable when the eigenvalues are negative and unstable 
when they are positive. The form of the trajectories near to the singular point is shown in the simulation results of Figure 
2.1, for the transfer function 1/(s* + 2.5s + 1) which has eigenvalues of -0.5 and -2, with corresponding eigenvectors 
of (1,-2)" and (1,-0.5)". 


Arrows are typically placed on the trajectories showing the direction of motion which will be towards the singular point for 
a stable node, as shown in the figure, and away from it for an unstable one. Trajectories starting on an eigenvector remain 
on it as is clearly seen in the figure. For the stable node all trajectories, apart from that starting on the first eigenvector tend 
towards the node along the second eigenvector, i.e the one with the smallest slope. The plot showing trajectory motions 
from different initial conditions is known as a state or phase portrait, although the latter name is often reserved for the 


special case where x; = x2, that is a = 0 and b = 1 in equation (2.3) 


Eigenvector 


x2 


Figure 2.1 Phase Portrait for a Node 


The singular point in case (ii) is a saddle point which is only theoretically reachable from initial conditions on the eigenvector 
of the stable solution associated with the negative eigenvalue. Any small perturbation will cause the trajectories to diverge 


as shown by the phase portrait of Figure 2.2 and move outwards as shown by the arrows. 
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Phase Portrait for Saddle Point 


Figure 2.2 Phase Portrait for a Saddle Point 


A phase portrait for case (iii) is shown in Figure 2.3 where the singular point is a focus. The focus is stable when the real 
part is negative, so that trajectories spiral towards it as marked in the figure, and unstable when the real part is positive. 


Four trajectories are shown in the figure. 


Phase Portrait for Focus 


Figure 2.3 Phase Portrait for a Focus Showing Four Trajectories 


In the final case (iv) the singular point is a centre, with the oscillatory motion producing concentric trajectories, as shown 
in Figure 2.4. Their size depends on the initial conditions, which in a physical situation, such as an ideal oscillating 


pendulum, is defined by the initial input energy. 
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Phase Portrait for Centre 


Figure 2.4 Phase Portrait for a Centre 


2.2.2. The Nonlinear Case 


In general it may not be possible to obtain analytical solutions to even second order nonlinear differential equations so 
the value of the phase plane approach in early work was to allow approximate solutions to be obtained using graphical 
techniques for sketching phase plane trajectories. The basic approach used was to determine the slope at a sufficient 
number of points in the state plane to allow a picture of the motion to be obtained starting from any initial conditions in 


the x, - x, state plane. From equation (2.2) a trajectory will have a slope r when it crosses the curve 


rP(x1,X2) + O(M1,xX2) = 0 (2.5) 
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By selecting a range of values of r, drawing the corresponding curves and marking the slope r with a short arrow at which 
a trajectory crosses, allows phase portraits to be sketched. This is usually known as the method of isoclines. Sometimes 


the lengths of the arrows are drawn in proportion to the velocity at the point and the plot is then called a vector field plot. 


Typically a second order nonlinear differential equation representing a control system with smooth nonlinearity can be 


written as 
x + f(x,x) = 0 


and if this is rearranged as two first order equations, choosing the phase variables, that is the output and its derivative, as 

the state variables, then x, = x, x. = x and the equation can be written as 
xX = X2 (2.6) 
X2 = — f(x1,x2) 

which is a special case of equation (2.2). It may have multiple singular points but they will all lie on the x, axis, that is 


for x, = 0. 


A phase portrait can therefore be obtained for any function, f, using the method of isoclines, or possibly other methods 
which can be found in the literature [2.1]. Digital simulation methods now enable phase plane trajectories to be easily 
displayed and they can often provide a clearer picture of the system behaviour than time response plots for x, and x,. 
Many situations in science and engineering may be modeled approximately by equation (2.1) or the special case equation 
(2.6). Some good examples can be found in the recent book given in reference [2.2] whilst here just one well known 


example is taken. 


2.2.3. An Example 


Many investigations using the phase plane technique were concerned with the possibility ofa nonlinear differential equation 
having a limit cycle solution, the strict name for an oscillation in a nonlinear system. When a limit cycle exists it results 
in a closed trajectory in the phase plane and typical of such investigations was the work of Van der Pol on oscillators. He 


considered the equation 
x-w—-x)xtx=0 (2.7) 


where // is a positive constant. The phase plane form of this equation can be written as 


X1 = Xo 

x2 =— f(m,x) = wd — xi) X2 — x1 

The slope of a trajectory in the phase plane is 
dx, x Wd — x1) x, — Xx, 


dx, x x) (2.8) 
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It has one singular point at the origin (0, 0) about which the linearised equation is x — x + x = 0, so that the singular 


point depends upon the value of 1. It is an unstable focus for y < 2 and an unstable node for yu > 2. 
All phase plane trajectories have a slope of r when they intersect the curve 
rx. = wd - xt) x2 — x (2.9) 


Figure 2.5 shows phase portraits obtained using isocline sketching and Figure 2.6 shows simulation results obtained from 
the initial conditions of (-3, 3) and (-0.1, 0.1) for u = 0.1 and 1. It can be seen, as expected, that the resulting limit cycles 
are dependent on u with that for the small value of 1 being almost elliptical. They are reached from outside for the first 
initial condition and from inside for the second initial condition. A limit cycle has this feature that it is reached from 
different initial conditions and in this case it is stable as trajectories tend to it from both without and within. In contrast 
if u = 0 then one has an oscillation, not a limit cycle, and the magnitude of the oscillation is determined by the initial 


conditions. Since the system has no damping there is no loss of energy and its value is set by the initial conditions. 


(d)u=5.0 


Figure 2.5 Phase portraits of the Van der Pol equation from isocline sketching 
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Simulation results for Yan der Pol equation with mu= 0.1 and 1.0 


Figure 2.6 Phase plane plots of the Van der Pol equation for two values of 


2.3. The Phase Plane for Systems with Linear Segmented Nonlinearities 


It is an advantage when using the phase plane approach to have nonlinearities which are described by linear segmented 
characteristics. This is because it results in a phase plane which can be divided up into different regions with different linear 
differential equations describing the motion in each region. Although the differential equations are linear the mathematical 
phase plane equations are only simple enough for calculating analytical solutions by ‘hand’ in a few cases. To illustrate 
these basic concepts three examples are considered in the next section. Simple dynamics are used in some cases to allow 


analytical solutions for the motion but also shown are some simulations with other dynamics. 
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2.3.1 Example 1 —- Nonlinear Output Derivative Feedback 


The Simulink diagram for this example is given in Figure 2.7, where the only nonlinearity is in the derivative of output 
feedback path. To plot the phase plane using Simulink the X-Y oscilloscope can be used or alternatively the output and 
its derivative can be fed to simout blocks for plotting after the simulation run has finished. The simout blocks must be 
put in the array mode and the time vector is available in tout. Often using the default integration algorithm in Simulink 
waveforms do not appear smooth so one may have to use an algorithm with a small fixed step size or limit the step size 
in an algorithm. Results from the X-Y scope can be read to the display accuracy. More accurate results can be obtained 
from the simout blocks but plotting the data may not be straightforward as one may have unequal vector lengths of data 


stored when a default integration algorithm is not used. 
If the nonlinearity is replaced by a gain K then the characteristic equation is 


s+(K-06)s+1=0 (2.10) 


To Wornkspacet 


Figure 2.7 Simulink diagram for system with dead zone feedback 


It is clear that for K < 0.6 the linear system is unstable and stable for K > 0.6. For the dead zone the gain is zero for a 
small sinusoidal signal but approaches unity for a large sinusoidal signal. Thus, one expects the effect of the dead zone 
nonlinearity to change the damping in some way such that the system is unstable for small signals and stable for large 
signals. This suggests the existence of a stable limit cycle. Figure 2.8 shows phase plane plots starting from (-1, 0) for the 


linear second order system with characteristic equation 
s+2f54+1=0 (2.11) 


for various values of €, the damping ratio. For the unstable case with ¢ = -0.05 the trajectory slowly spirals outwards 
and for the stable responses they converge to the origin in a spiral manner for & < 1, when the origin is a focus and 
directly with no overshoot for € = 1, when the origin is a node. The mathematical expressions for these curves are quite 
complicated, so to sketch them on a phase plane the method of isoclines or some special methods may be used [1]. The 
dead zone is set at levels of +0.5, so that the trajectories describing the motion are those of a linear second order system 
with a damping ratio of -0.3 for |x,| < 0.5 and 0.2 for |x,| > 0.5. Phase plane plots have been obtained by simulation and are 


shown from initial conditions of (-8, 0) and (0.1, 0) in Figure 2.9. The resulting limit cycle is clearly shown in the figure. 
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Damping ratio zeta 


Figure 2.9 Phase plane plots for system of Figure 2.7 


2.3.2. Example 2 — Relay Position Control 


Consider a basic relay position control system with nonlinear velocity feedback having the Simulink diagram shown in 
Figure 2.10. The two relays in parallel simulate a relay with dead zone and hysteresis which has the characteristic shown in 
Figure 2.11, where for positive inputs the switch on level is 6 + A and the switch off level is 0 — A. The relay parameters 


are set to give an output, h, of + 1, and 0. 
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Saint Integrator 


To Workspace 


Saturation 


Figure 2.10 Simulink diagram of the relay control system for example 2 


Figure 2.11 Relay with dead zone and hysteresis 


It is assumed initially that the hysteresis in the relay is negligible (ie. A = 0) and that there is no saturation in the velocity 

feedback path.. Denoting the system position output by x; and its derivative x; by x2 then the feedback to the relay input 

is assumed to be — x: — Ax2. The relay output of + 1 or 0 is equal to x: /K, where in the Simulink diagram both K and A 

are chosen equal to 0.5. Taking the dead zone of the relay +0 to be equal to +1, the motion of the system is described by 
K if=a=Aee> 1 

x=90 ifl-xm—Ami< 1 (2.12) 
=Kit=m = Ae <1 

Thus the equation of motion in the phase plane changes at the lines x, + Ax. =+ 1 and these can be drawn on the plot, 

as shown in Figure 2.12, to divide up the phase plane into the three regions where the motion is described by the above 


three simple linear second-order differential equations. 
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-2 


-3 -1 
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The dotted lines show the 
change in the switching 
boundaries due to 
saturation 


Figure 2.12 Switching lines in the phase plane 
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is obtained by writing it in the form 


a (2.14) 
dx, 


and integrating with the initial conditions for x, and x, of x,,and x 


49 Lespectively, to give 


x2 — X30 = 2K(m — x0) (2.15) 


where x10 and X29 are the initial values of x, and x2. Since equation (2.15) describes a parabola, which for the special 
case of K = 0 has the solution x2 = x29, it is easy to calculate the system’s response from any initial condition (x0, x20) 


in the phase plane. 


Figure 2.13 shows the response from (-2, 0) with A = K = 0.5 as given in Figure 2.7. The initial parabola, with K = 0.5, 
starts from A and meets the first switching boundary x: + 0.5x2 =— 1 at B; the ensuing motion is horizontal, that is, 
at constant velocity, until the second switching boundary x; + 0.5x2 = 1 is reached at C. The ensuing parabola, with 
K=-0.5, meets the same switching boundary again at D. The next motion is again at constant velocity until the point E is 
reached on the first switching boundary. At this point the following parabolic motion, with K = 0.5, brings the trajectory 
straight back to the switching boundary. The situation has therefore arisen where the motion at either side of the switching 
line is directed back to it. Thus, the resulting motion, which is known as a sliding mode, is directed along the switching 
line to come to rest at (-1, 0). In theory switching takes place at an infinite rate but in practice the relay will have some 


hysteresis, which will decrease the switching rate. 


By solving the relevant equations it can be shown that B = (-1.39, 0.781), C = (0.610, 0.781), D = (1.14, -0.281) and 
E = (-0.86,-0.281). If the first integrator is replaced by the transfer function 1/(s + a) to provide some damping then the 
analytical solution is not as easy. The motion in the phase plane when the relay output is zero is still linear, but now with 
a slope of -1/a, and when the relay output is +1 it is no longer parabolic. Simulation responses are also shown in Figure 
2.13 for values of a equal to 0.1 and 0.3. Responses from any other initial conditions are obviously easy to find, but, from 
the responses shown, several aspects of the system’s behavior are readily apparent. In particular the system is seen to be 
stable since all responses will move inward, possibly with several overshoots and undershoots, and will finally slide down 


a switching boundary to +1. Thus a steady-state error of unit magnitude will result from any motion. 


Finally a word of warning that difficulties may arise in simulating systems with sliding, since if a variable step length 
integration algorithm is used it may reduce the step length to zero and stop in the sliding motion. Many simulation 
languages have special routines for taking care of the discontinuous sliding motion, but good results can usually be obtained 


by using a fixed step length integration algorithm or by incorporating a small hysteresis in the relay. 
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Figure 2.13 Initial condition response for system with 


When the velocity feedback signal saturates, that is, when |Ax2|> /, the input signal to the relay is — x; +h. Thus, if h is 
chosen equal to 0.5 when |x, |> 1 the switching boundaries become vertical lines as shown dotted in Figure 2.12, although 
the equations describing the motion between the boundaries remain unaltered. Therefore for a large step input of 1, the 
equivalent of starting from the initial condition (-r, 0), the response will become more oscillatory when the velocity 
saturates as illustrated in Figure 2.14 for values of r equal to 6 and 3 and saturation of the velocity feedback signal at +0.5 
for inputs greater than +0.5. The parabolic curves are clearly seen, since the damping has been taken as zero, as also is 


the fact that the switching lines have become vertical for |x2|> 1. 


Figure 2.14 Phase plane plots for step inputs of 6 and 3 with velocity saturation but no hysteresis or damping. 


Again if there is no velocity saturation in the feedback path but the relay has a finite hysteresis, A, then the switching 
boundaries are again straight lines of slope -1/A but different according to whether x, is positive or negative.In the general 


case they are 


x +Ax.=-O+A and x +Ax. = 6+ A for x2 positive, and (2.16) 


x +Ax=-dO-A and x +Axm = 0 —A for x2 negative. 
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It is easily shown from symmetry considerations, when the above lines are drawn on a phase plane, that a limit cycle 
exists for the system with no damping with a maximum value of x. = A/A. Figure 2.15 shows simulation results for 
0 = land A = A = 0.5. The responses from -1.51 and -4 are shown converging to the limit cycle from within and 
without. If the system were stabilized by replacing the first integrator by the transfer function 1/(s + a) then the motion 


could come to rest, dependent on the input step magnitude, with the relay input in the range +1.5. 


Figure 2.15 Initial condition responses resulting in a limit cycle 
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2.3.3. Example 3 — Position Control with Torque Saturation 


Figure 2.16 shows a Simulink diagram for a position control system with no viscous damping and with nonlinear effects 


due to torque saturation and Coulomb friction. 


Gain Saturation 


To Wornspacet 


Figure 2.16 Simulink diagram of control system for Example 3 
The differential equation of motion in phase variable form has 
x2 = f(— 1) — 0.5sgn (x2) (2.17) 
where jf; denotes the saturation nonlinearity and sgn the signum function, which is +1 for x. > 0 and -1 for x < 0. 


There are six linear differential equations describing the motion in different regions of the phase plane. For x2 positive, 


equation (2.17) can be written 


x + f(x) + 1/2 = 0 (2.18) 
so that for 
(a) x2 + ve,x < — 2, one has x, = x2,xX2 = 3/2, a parabola in the phase plane. 
(b) x2 + ve, — 2 <m S< 2, one has x) = x2,%2 +x + 1/2 = 0, a circle in the phase plane. 


(c) x2 + ve,x, > 2, one has x, = %2,x%. =— 5/2,a parabola in the phase plane. 


Similarly for x2 negative, 


(d) x2 — ve,x1 < — 2, one has x1 = x2,X2 = 5/2, a parabola in the phase plane. 
(e) Xx. — ve, — 2 <x < 2, one has x = x2,xX. + x; — 1/2 = 0, a circle in the phase plane. 
(f) x2 — ve,x; > 2, one has x, = x2,xX2. =— 3/2, a parabola in the phase plane. 
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Because all the phase plane trajectories are described by simple mathematical expressions, it is straightforward to calculate 


specific phase plane trajectories. 


Figure 2.17 Response from initial condition (-6, 0) for Example 3 


Careful examination of Figure 2.17, which shows the simulation result for a response from (-6, 0), reveals the changes in 


the trajectory shape when crossing the lines x, = +2. 


2.4 Conclusions 


This chapter has covered the basic concepts in the analysis of second order systems using the phase plane approach. 
Since many simple control systems can be approximated by nonlinear second order differential equations the approach is 
powerful. It can be used when more than one nonlinear element exists and is particularly useful when the nonlinearity can 
be approximated by linear segmented characteristics. Examples have been given using a double integrator plant transfer 
function with a constant input, which results in simple mathematical expressions for the phase plane trajectories, so that 


certain features can be easily illustrated. 


Mathematical expressions for the trajectories can always be obtained, since they are given by the solutions of linear 
differential equations when the nonlinearities are linear segmented characteristics, but they are often quite complicated. 
Phase plane plots can easily be obtained with modern simulation languages but the facility to sketch a phase portrait is 
an important aid to understanding the system behaviour. The results can also be useful for comparing results obtained by 


the approximate DF method, to be introduced in the next chapters, when it is applied to second order systems. 
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3 The Describing Function 


3.1 Introduction 


The describing function, which will be abbreviated DF, method was developed simultaneously in several countries during 
the 1940s. Engineers found that control systems which were being used in many applications, for example gun pointing 
and antenna control, could exhibit limit cycles under certain conditions rather than move to a static equilibrium. They 
realized this instability was due to nonlinearities, such as backlash in the gears of the control system, and they wished to 
obtain a design method which could ensure the resulting systems were free from limit cycle operation. They observed that 
when limit cycles occurred the waveforms at the system output were often approximately sinusoidal and this indicated to 
them a possible analytical approach, namely to assume that the signal at the input to the nonlinear element in the loop 
was a sinusoid. Since then there have been many developments in terms of both using the DF concept for other types of 
signals and the problems, or phenomena, which they can be used to study. More will be said on these aspects later but 
we begin by considering the initial problem of investigating the possibility of a limit cycle in a feedback system using the 


DF or S (sinusoidal) DF as it is often named. 


Consider the feedback system shown in Figure 3.1 containing a single static nonlinearity n(x) and linear dynamics given 
by the transfer function G(s) = G.(s)Gi(s). Ifa limit cycle exists in the autonomous system, that is with r(t) = 0, with the 
output c(t) approximately sinusoidal, then the input x(t) to the nonlinearity might also be expected to be near sinusoidal. 
If this assumption is made the fundamental output of the nonlinearity can be calculated and conditions for the sinusoidal 


self-oscillation found, if the higher harmonics generated at the nonlinearity output are neglected. 


This is the concept of harmonic balance, in this case balancing the first harmonic only. This approach had previously 
been used by physicists to investigate such aspects as the generation of oscillations in electronic circuits; the Van der Pol 
oscillator mentioned in the previous chapter being an example. The DF of a nonlinearity is therefore defined as its gain 
to a sinusoid, that is the ratio of the fundamental of the output to the amplitude of the sinusoidal input. Since the output 


fundamental may not be in phase with the sinusoidal input the DF may be complex. 


“ote H Heh 


4 


Figure 3.1 A simple nonlinear feedback system 


3.2 The Sinusoidal Describing Function 


Assume that in Figure 3.1 x(t) = a cos 0, where 9 = wt and n(x), where and n(x) is a symmetrical odd nonlinearity, 


then the output y(t) will be given by the Fourier series. 
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CO 
y(@) = >¥ a,cosnO + b, sin nO (3.1) 
n=0 
where d, = b, = 0 for n even, and in particular 
27 
a = (1x) [" yO) cos 0 ad (3.2) 
and 
20 
bi = (ia) [" yO) sin 0 a9 (3.3) 
The fundamental output from the nonlinearity is a, cos 9 + b; sin @, so that the DF is given by 
N(a) = (a — jhi)/a (3.4) 
which may be written 
N(a) = N,(a) + jNq(a) (3.5) 
where 


N,(a) = a/aand N,(a) =— b\/a (3.6) 


Alternatively, in polar co-ordinates 


N(a) = M(ae™ (3.7) 
where 

M(a) = (ai + bi)" la (3.8) 
and 

v(a) =— tan '(b)/a). (3.9) 


If n(x) is single valued it is easily shown that b; = 0 and from symmetry considerations as n(x) is odd, that 
m2 
a =a) i (0) cos 040 (3.10) 
0 
giving 
a2 
N(a) = (aila) = (4/az) ‘i y(0) cos 040 (3.11) 
0 
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Although equations (3.2) and (3.3) are an obvious approach to the evaluation of the fundamental output of a nonlinearity, 
they are somewhat indirect, in that one must first determine the output waveform y(@) from the known nonlinear 
characteristic and sinusoidal input waveform. This is avoided if the substitution 6 = cos | (x/ a) is made; in which case, 


after some simple manipulations, it can be shown that 


a, = (2/a) J smpCp ax (3.12) 
bi = (2/az) f weteldk (3.13) 


The function p(x) is the amplitude probability density function of the input sinusoidal signal and is given by 

p(x) = (May(a@® — x’ (3.14) 
and the nonlinear characteristics np(x) and n,(x), called the in-phase and quadrature nonlinearities, are defined by 

Np(x) = [m(x) + m(x)]/2 (3.15) 
and 

Ng (x) = [N2(x) — m(x)]/2 (3.16) 


where 1(x) and n2(x) are the portions of a double valued characteristic traversed by the input for x > 0 and x < 0 


respectively. 


For a single-valued characteristic, m(x) = m2(x), so that n)(x) = n(x) and n,(x) = 0 and integrating equation (3.12) 


by parts gives 
a, = (4/)n(0*) + (alan) |" n'(x)(@ — x?) dx (3.17) 
0 
where n' (x) = dn(x)/dx and n(0°) = lim n(é);a useful expression for obtaining DFs for linear segmented characteristics. 


3.3. Some Properties of the DF 


An advantage of the formulation of equations (3.12) and (3.13) is that they easily yield proofs of some interesting properties 


of the DF for symmetrical odd nonlinearities. These include the following: 


1. For a double-valued nonlinearity the quadrature component N, (a) is proportional to the area of the 
nonlinearity loop, that is: 


N,(a) =— C/ a’) (area of nonlinearity loop) 


2. For two single-valued nonlinearities ng(x) and ng(x), with ne (x) < ng(x) for all 0 < x < b, then 


N2(a) < Ng(a) for input amplitudes a less than b. 
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3. For the sector bounded single-valued nonlinearity that is kjx < n(x) < k2(x) for all 0 < x < b then 
k, < N(a) < kp for input amplitudes a less than b. This is the sector property of the DF and it also applies for 
a double-valued nonlinearity if N(a) is replaced by M(a). 


When the nonlinearity is single valued, it also follows directly from the properties of Fourier series that the DF, N(a), 


may also be defined as: 


1. The variable gain, K, having the same sinusoidal input as the nonlinearity, which minimizes the mean 


squared value of the error between the output from the nonlinearity and that from the variable gain. 


2. The covariance of the input sinusoid and the nonlinearity output divided by the variance of the input. 


It can also be shown that the gain K is the best transfer function model in the mean squared error sense. That is the 


optimum G(s) is K. 
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3.4 The Evaluation of some DFs 


To illustrate the evaluation of the DF for some specific nonlinearities a few simple examples are considered below. 


3.4.1 Cubic Nonlinearity 
For this nonlinearity n(x) = x and using equation (3.10) one has, 


a= (ain) [~ (acos 6)* cos 6d 
0 


n/2 
= (Alaa f *6d0 
a (3.18) 
_ 3 (77/3 , cos26 , cos4é sas & 
= (4/n)a f & + 5 + 9 )ao = 3a /4 giving 
Na) = 3a°/4 


Alternatively from equation (3.12) and using the fact that n(x) is odd, one has 
a= (41a) | °x‘p(aax 
0 


The integral 4, = ‘i x" p(x)dx is known as the n* moment of the probability density function and for the sinusoidal 


distribution with p(x) = (1/z)(a’ — x°)"”, fm, has the value 


0 for n odd 
Ln _ le = = 
a (n— 1) @—3) iiss ] for n even 
n (n-2) 2 
Therefore N(a) = (Ala) a = 3a°/4 as before. 


3.4.2 Saturation Nonlinearity 


The DF can also be found by taking the nonlinearity input as asin 0, in which case for the ideal saturation characteristic 
shown in Figure 3.2, with unit slope in the linear regime and saturation for an input above 0 > 0.5, the nonlinearity 


output waveform y(@) is as shown in the same figure. 


Ideal saturation nonlinearity 1 


output 


Saturation level 0.5 


L L A 1. . A 1 L 1 n 1 1 L 1 1 
cP ae) Sas Ste ales Q a2 04 O06 O8 1 it) 1 2 3 4 6 6 i 8 9 10 


Figure 3.2 Saturation nonlinearity, input sinusoidal and output waveforms 
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Again, because of the symmetry of the nonlinearity the fundamental of the output can be evaluated from the integral 


over a quarter period, so that for a linear regime slope of m and saturation at 6 gives 
4 m2 
N(a) = _ f'y(@)sin 6a9 
ar Jo 
which for a > 0 gives 
4 a 2 al2 . 
N(a) = — masin 0d@ + | mdsin6d0 
art 0 a 
with @ = sin '6/a. Evaluation of the integrals gives 
N(a) = (4min)| © sBm 4! 2 pos @| 
2 4 
which on substituting for 0 gives the result 
N(a) = (m/z)(2a@ + sin2a). 
Since for a < 6 the characteristic is linear giving N(a) = m, the DF for ideal saturation is mN;(6/a) where 


mN,(6/a) = \ fora <d (3.19) 


(1/m)[2a + sin2a] for a> 


Alternatively one can evaluate, N(a), from equation (3.17), which yields, as the nonlinearity slope is zero for inputs greater 
than 6 


8 
N(a) = ala = (4la’ x) f m(a? =) ae 
0 
Using the substitution x = asin J, this gives 
oe (amix) f° cos’ 0d@ = (m/x)(2a + sin 2a) as before. 
0 


3.4.3. Relay with Dead Zone and Hysteresis 


The characteristic of a relay with dead zone and hysteresis, for 0 = 1 and A = 0.5, is shown in Figure 3.3 together with 


the corresponding input, assumed equal to acos @, and output waveforms. 
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Figure 3.3 Relay with dead zone and hysteresis, and input and output waveforms 


Using equations (3.2) and (3.3) over the interval — 7/2 to z/2 and assuming that the input amplitude a is greater than 


0 + A, gives 
B 
a = (21) | heos 6d0 = (2h/z)(sinB + sina) 


where @ = cos ‘(6 — A)/a and B = cos '(6 + A)/a, and 


B — 
bh = (In) f hsin 6d@ = nin(O* A) _ (0 = 


) = 4hA/az. Thus 
a 


Nia = 2 f[a? — (6 + Ay] + [a - 6 - AY} - BE 6.20) 


2 
ag 


Do you want to 
make a difference? 


Join the IT company that 
works nard to make life 


easier. 


www.tieto.fi/careers 


Knowledge. Passion. Results. 


Download free ebooks at bookboon.com 


An Introduction to Nonlinearity in Control Systems The Describing Function 


For the alternative approach one must first obtain the in-phase and quadrature nonlinearities which are shown in Figure 


3.4. Using equations (3.12) and (3.13) one obtains 
é+A a 
a = (41a) [ x(hI2)p(x) dx + J xhp oa 
6-A O+A 
_ an @ =(04 Ay}? di [a - (6 = Ay FP} 
ar 
6+A 
b = (Alar) [ (hi2) dx = 4hA/az = (Area of nonlinearity loop) /az. 
d—A 


The DFs for other relay characteristics can easily be found from this result or worked out directly. The latter is very 
straightforward when the relay has no dead zone as its output is a square wave of amplitude th, which has a fundamental 
component of amplitude 4h/7, either in phase with the input for A = 0, or lagging when A is finite. Using equation 3.20 


for no hysteresis, A = 0; for no dead zone, 6 = 0; and for an ideal relay, A = 0 = 0. 


o/p o/p 
1.0 
- “AU 
05 15 Wp 0.5 
In phase characterisitic Quadrature characteristic 


Figure 3.4 In phase and quadrature nonlinearities 


It is possible to model a typical double valued nonlinearity like the above as a nonlinear function n(x,x) and the describing 
function theory based on this can be found in reference 3.1. Unfortunately this description is not correct for any input 
signal, x, and it has been used to obtain erroneous results for x a random signal. Two tables of DFs for a sinusoidal 
input are given in the Appendix 3.11, one for single valued nonlinearities and the other for double valued nonlinearities. 


Extensive use is made in the tables of the DF for saturation, N,(0/a), defined in equation (3.19) 


3.5 Nonlinear Models and DFs 


Nonlinear models, which consist of linear segments with multiple break points, can often be modeled from the nonlinearities 
in (b) and (c) above and linear gains, usually connected in parallel. For example, the dead zone characteristic which for a 
positive input has a zero output until a value of 6 and then a slope of unity is easily shown to be modeled by a unit gain 
in parallel with a saturation characteristic with slope m = -1, as illustrated in Figure 3.5. Also a quantized characteristic 


can be modeled from relay characteristics with dead zone and no hysteresis in parallel. 


Download free ebooks at bookboon.com 


44 


An Introduction to Nonlinearity in Control Systems The Describing Function 


Figure 3.5 Modelling of an ideal dead zone characteristic from a unit gain and ideal saturation 


Since it is easily shown that the DF of two nonlinearities in parallel is equal to the sum of their individual DFs, the 
DFs of linear segmented characteristics with multiple break points can easily be written down from the DFs of simpler 
characteristics. Several procedures are available for obtaining approximations for the DF of a given nonlinearity either by 
numerical integration or by evaluation of the DF of an approximating nonlinear characteristic defined, for example, by a 


quantized characteristic, linear segmented characteristic or Fourier series. 


3.6 Harmonic Outputs 


It is sometimes of value to know the magnitudes of the other harmonic outputs from a nonlinearity with a sinusoidal 


input. These can be obtained from the usual Fourier series results for equation (3.1) of 
20 
an = (1/7) f n(x) cos n@d@ (3.21) 
0 
and 
20 
b, = (1/z) [rc sin n6d0 (3.22) 
0 


In terms of the amplitude probability density function for a sinusoid it can be easily shown that these expressions become 


fm = 2 f n(x) T.(xla)p (x) de (3.23) 
and 
b, = 2 f n@Da(ala) paar (3.24) 


where 7,,(x/a) = cos|n cos | (x/a) | and D,(x/a) = sin|n cos | (x/a)| 
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Since for single valued nonlinearities b, = 0 the more important result is for a,. Here T (x/a) is the Chebyshev polynomial 


of order n and its value for any n can easily be found from the recursion relationship 
Tr+ (x) = 2xTh(x) — Th-1 (x) (3.25) 
with 7)(x) = land 7i(x) =x 


3.7. Sine plus Bias DF and the IDF 


So far only DFs for nonlinearities with odd symmetry have been considered. For many control engineering situations 
this may be satisfactory, as for example in the consideration of torque saturation in a motor. Many situations do occur, 
however, where nonlinearities do not have odd symmetry and also even when they do the operation may not be about 
the odd symmetric axis. Thus, to handle the general situation for a single valued nonlinearity with characteristic n(x) 
one needs to consider an input consisting of a sinusoid plus bias, that is y + acos @ and evaluate both the sinusoidal and 
bias output of the nonlinearity. The describing function for the nonlinearity, known as the SBDF, has two components 
N,(a, 7) and Na(a, y) the gains to the bias and sinusoid, respectively. It is easily shown following the procedures of section 
3.2 that 


Nj(a,7) = (aol7) = (Uy) f not Ypewrax (3.26) 
and 
Nala, v) = (a:/a) = (2/a) fre + ¥)(x/a) p(x) dx (3.27) 


Thus for the input y + acos @ the output is yN;(y,a) + Na(y,a)acos @. The important point about the SBDF model is the 
interaction between the signals caused by the nonlinearity, so that the gain to both signals is affected by a change in one 
of them. If a is constant and y varied the bias output of equation (3.26) can be plotted against y to give a new nonlinear 
characteristic. This is often called a modified nonlinearity and denoted by n(a, y) rather than ao. The slope of the modified 


nonlinearity at any input value y gives the gain g(a, 7) to any small unrelated signal in the presence of a and y. That is 


g(a,y) = a fn + y)p(x)dx (3.28) 


In particular, the gain to a small unrelated signal in the presence of x alone is g(a,0), usually called the incremental 


describing function (IDF) and denoted by Nj,(a), and is given by 


a 


Ny(a) = g(a,0) = f n' (x)p (x) dx (3.29) 


where n'(x) = dn(x)/dx. Equation (3.29) is the bias output from the nonlinearity n' (x) with input x. For any single valued 


nonlinearity with n(0*) = 0, equation (3.17) gives 


aN) OM i. n' (xe) (a? — «de 
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which on differentiating both sides with respect to a gives 
2aN(a) + a’dN(a)/da = 2 fan cop ends 
The left hand side is 2aNi,(a), so that 
Niy(a) = N(a) + (a/2)N' (a) (3.30) 
where N'(a) = dN(a)/da. 


3.8 Conclusions 


This chapter has covered the definition, some properties, and approaches to evaluating the DF of a nonlinearity for a 
sinusoidal input. The approach has then been extended to cover the situation of a sinusoidal plus bias input. The next 


chapter will show how these results can be used in the approximate analysis of a simple nonlinear feedback system. 
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3.11 Appendix -Tables of Describing Functions 


Two tables of DFs, or SDFs, are given below, one for single valued nonlinearities and the other for double valued 
nonlinearities. Extensive use is made in the tables of the describing function for saturation, N;(0/a), defined in eqn. 
(3.19). More comprehensive tables are given in other books [3.1, 3.2] for other describing functions, such as the SBDF 


and IDF of section 3.7, and for other inputs. 


Download free ebooks at bookboon.com 


47 


Nonlinearity 


General quantizer 


An Introduction to Nonlinearity in Control Systems 


Comments 


a<6, 


The Describing Function 


N,(a) = a,/a 


N,=0 


M 
he bu+i>a> dy N, = (4/a?x) Y h,,(a? — 53)” 
A. hs ' m= 
Ar | 
8, 8,53 5, 
Uniform quantizer a<6é N, =0 
hy =h,=+++h M 
Sm = (2m — 1)5/2 (2M + 1)6 >a>(2M — 1)6 N, = (4h/a?x) > (a? — 262)? 
Ban 
n=(2m-1)/2 
Relay with dead zone a<6 N,=0 
a a>o N, = 4h(a? — 5°)"?/a*x 
r) 
Ideal relay N, = 4h/ax 
h 
Preload N, = (4h/ax) + m 
A m 
General piecewise linear a<6é, N, = (4h/ax) + m, 


Suri > «> dy Ny = (4h/am) + my 


M 
+ 2 (m; — mj,,)N,(5,/a) 


Ideal saturation N, = mN,(6/a) 
ge = 
8 
Dead zone N, = m[1 — N,(6/a)] 
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a<éd 


a>éd 


a<6 


a>é 


a<6 


a>o 


Limited field of view 


a<é 
a a>é 
) 


y=x m > —2 Tis the gamma function 

y=x" n odd integer. y,, is the nth 
moment of the amplitude 
distribution function of a 
sinusoid 

y=x" n even integer 


Harmonic nonlinearity J,(ma) is the Bessel function 
y =A, sin mx of order 1 


y =A, sinh mx (ma) is the modified 
Bessel function of order 1 
Exponential saturation 1,(ca) is the modified 
y=1-e*% Bessel function of order 1. 
S,(ca) is the modified 
Struve function of order 1 


cx K(k) and E(k) are complete 
Pr as tN elliptic integrals 
Brees u=ca/(1 +c? 


The Describing Function 


N, = —m,N,(6;/a) + (m, — m.)N,(42/a) 


+m, 


N,=0 


N, = 4h(a? — 8?)'?/a,0 + m — mN,(6/a) 


N, =m, 


N, = (m, = m,)N,(6/a) +m, 
+ 4h(a’? — 5°)'?/a?x 


N, = 4h/ax 
N, = 4h/[a — (a? — 5°)'")/a?x 


N, = (m, + my)N,(8/a) 
—mN[(m, + m)6/m,a] 


N, =m, 


N, = m,N,(6/a) — 4m,6(a? — 5?)"?/a*n 


C(m +a"! 


»~ 2=IT[(3 + m)/2\r[(1 + my/2] 
2 T[(m + 2)/2]a""! 


“Vx P[(m + 3)/2] 
_ n(n — 2)(n — 4)-+-3 : 
P (n+ (n= 14 


2 
= i tat 


N 


N 


4 n(n-2)---2 


ge Es So 


N, = 2AmJ ;(ma)/a 
N, = 2A m1 \(ma)/a 


N,= (2/a)[J (ca) — S,(ca)] 


N, = (4/m)[(—u/c?a*)K(u) + (1/au)E(u)] 


Table 3.1 Describing Functions of Single Valued Nonlinearities 
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Nonlinearity Comments N,(a) = a,/a Na) = —b,/a 
On-off relay with a<A N,=0 N,=0 
hysteresis 


¥ a>A N, = 4h(a? — A*)'?/a?x N, = ~4hA/a'x 


Friction-controlled 


backlash a>b/2 N, = {1 + N,[(a — 5)/a)}/2 N, = b(b — 2a)a*x 
Al 
5/2 
>A +(h/m) N oy, (24) +m(224)] N, = —4hA/a® 
h a m p=alN, iq q=7 ax 
es 2 ma ma 
4 
h+mA —mA 
/ PL o>asm n= Fl 2—n,(*2tA) 4 (MY) nanan 
4 
i] m,>m, 
h Am, * a jee 
im Om —m) n,-( 2 N, ma Ny = —4hA/a?x 
we Am, h—ma 
(m, —m)) sey 14 bles 
2 am 
a>A N, =m + 4h(a? — A?)!?/a?x y= ~4hA/a'n 
4 
On-off relay with a<d+A N,=0 N,=0 
dead zone and 2h ‘ sik 
hysteresis a>d+A Np, = an l2 —(6+4)7)" N, = —4hA/a’x 


ny 28 + [a - (6 - )?}"7} 
it) @ 


Table 3.2 Describing Functions of Double Valued Nonlinearities 
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4 Stability and Limit Cycles using the 
DF 


4.1 Introduction 


The problem of determining whether a linear feedback system is stable is simply one of determination of the roots of 
its characteristic equation. This can be done in several ways with a particularly elegant and convenient engineering one, 
since it can be done from knowledge of the open loop frequency response locus, being due to Nyquist. On the other 
hand the problem of determining precisely the stability of an autonomous nonlinear feedback system, even one simply 
containing one linear dynamic element and one nonlinearity, has not been solved, even though it has exercised the minds 
of mathematicians and engineers for many years. The DF provides an approximate method for answering the problem and 
will be addressed in this chapter, whilst some exact approaches, known as absolute stability methods, will be considered 
in chapter 8. The discussions will be restricted to the feedback loop structure of Figure 3.1. Since our initial investigations 
are concerned with the stability of the autonomous system, that is r(t) = 0, the two linear dynamic blocks are in series 
and can be represented by the single transfer function G(s). Further the position of the single static nonlinearity although 


assumed in the forward path could equally well be in the feedback path. 
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An early contribution to the problem was a conjecture by Aizermann. He conjectured that ifa symmetrical odd nonlinearity 
was confined within a sector defined by straight lines of slope k, and k, then any nonlinear system with a nonlinearity lying 
entirely within the sector would be stable provided the linear system was stable for gains between. k, and k,. Unfortunately 
provided the linear system was stable for gains between. k, and k,’ Thus, since use of the DF provides an approximate 


method for the determination of limit cycles, it also provides an approximate stability test for our simple feedback loop. 


4.2 Limit Cycle Evaluation 


To study the possibility of limit cycles in the autonomous closed loop system of Figure 3.1, the nonlinearity n(x) is replaced 


by its DF N(a). Thus, the open loop gain to a sinusoid is N(a)G(jw) and a limit cycle will exist if 
N(a)G(ja) =- 1 (4.1) 


where G(jw) = Gc(jw) Gi (jw). This condition means that the first harmonic is balanced around the closed loop assuming 
its passage through the nonlinearity is accurately described by N(a). Since G(jw) is a complex function of w and N(a) 
may be a complex function of a, a solution to equation (4.1) will yield both the frequency w and amplitude a of an 


assumed sinusoidal limit cycle. 


Various approaches can be used to examine equation (4.1) with the choice affected to some extent by the problem. For 
example, relevant factors might be whether the nonlinearity is single or double valued or whether G(jw) is available from 
a transfer function G(s) or as measured frequency response data. Typically the functions G(jw) and N(q) are plotted 
separately on Bode, Nyquist, or Nichols diagrams. Alternatively, stability criteria such as the Hurwitz-Routh or root locus 


plots may be used for the characteristic equation 
1+ N(a)G(s) = 0 (4.2) 
although here it should be remembered that the equation is appropriate only for s ~ jw. 


Figure 4.1 illustrates the procedure on a Nyquist diagram, where the G(jw) and C(a) =— 1/N(a) loci are plotted and 
shown intersecting at P for a = ad) and @ = Wo. The DF method therefore indicates that the system has a limit cycle with 
the input sinusoid to the nonlinearity, x, equal to a) sin(wot + 6), where & depends on the initial conditions and time 
origin. In practice the limit cycle will not be sinusoidal and ao is an approximation for the amplitude of the fundamental 
component of the limit cycle. Thus if one wishes to estimate the accuracy of the DF prediction for a limit cycle one should 
measure the amplitude of the fundamental not the peak amplitude of the waveform as is often done for convenience. 
When the G(jw) and C(a) loci do not intersect, the DF method predicts that no limit cycle will exist if the Nyquist 
stability criterion is satisfied for G(jw) with respect to any point on the C(a) locus. Obviously, if the nonlinearity has 
unit gain for small inputs, the point (— 1,j0) will lie on C(a) and it may then be used as the critical point, analogous 


to the situation for a linear system. 
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= 
Ca) = -1/N(a) 
“4h Gijw) 
Nyquist 


Figure 4.1 Limit cycle found from intersection of G(jw) and C(a) ona Nyquist plot 


When the analysis indicates the system is stable its relative stability may be indicated by evaluating its gain and phase 


margin. These can be found for every amplitude a on the C(a) locus, so it is usually appropriate to use the minimum values 


4.3 Stability of a Predicted Limit Cycle 


A further point of interest when a solution to equation (4.1) exists is whether the predicted limit cycle is stable. This is 
obviously important if the control system is designed to have a limit cycle operation, as in the case of an on-off temperature 
control system, but it may also be important in other systems. If, for example, an unstable limit cycle condition is reached 
the signal amplitudes will not become bounded but continue to grow. The stability of a limit cycle, provided only one 
solution is predicted, can in most cases be assessed by applying the Nyquist stability criterion to points on the C(a) locus 
at both sides of the solution point. If the stability criterion indicates instability (stability) for a point on C(a) with a < ao 
and stability (instability) for a point on C(a) with a > a,, then the limit cycle is stable (unstable). This result is known as 
the Loeb intersection criterion and in terms of the approximate DF analysis, is a necessary, but not a sufficient criterion 


for a stable limit cycle. 


To investigate the stability of a limit cycle one has to examine the perturbation equation. Consider the transfer function 


to be the ratio of two polynomials so that it can be written as 
G(s) = p(s)/q(s) (4.3) 
Then if D denotes the differential operator the differential equation for the autonomous system can be written 
q(D)x(t) + p(D)n(x(t)) = 0 (4.4) 
If the system has a limit cycle solution x * (¢) and this is perturbed by a small amount Ox (ft) then 


q(D)[ x * (t) + dx(O)] + p(D)n(& * (t) + dx(D) = 0 
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This can be written as 


II 
S 


q(D) [x * () + dx()] + pD)ne* (H) + PCD)n'(x * (Hdx(O) 
Thus on subtracting equation (4.4) one has the perturbation equation 
q(D) 0x (t) + p(D)n' (x * ()) dx) = 0 (4.5) 


This is a time varying equation for the perturbation, 0x(t), which must be stable for the limit cycle to be stable. It can be 
shown that a necessary condition for this equation to be stable is that the time averaged equation must be stable. Then if 
one uses the approximate describing function solution for x * (t), the time averaged value of the nonlinear term is the IDF 
of equation (3.29), so that a necessary condition for stability of the limit cycle is that the equation g(D) + p(D)N,,(a) = 0 


is stable. In terms of the system transfer function this is the characteristic equation 


1+ N,(a)G(s) = 0 (4.6) 
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4.4 DF Accuracy 


A specification for many control systems is that when perturbed by either a step input or disturbance they should return 
to equilibrium. For a nonlinear system this means checking that (a) the response will not go unbounded and (b) no limit 
cycle will occur. Condition (a) can be checked from linear methods using the sector bounds of the nonlinearity provided 
any input does not cause a bias at the nonlinearity input and (b) by the DF method. The DF method, however, can only 
do this approximately since it assumes any limit cycle will be sinusoidal at the input to the nonlinearity, which will never 
be quite true in practice. This means that the DF method could indicate incorrectly either the existence or non-existence 


of a limit cycle. It is therefore important to have some idea of the validity and accuracy of any DF result. 


If the DF method predicts a limit cycle then its validity can normally be checked by assuming this sinusoidal signal as the 
nonlinearity input and evaluating the signal fed back to the nonlinearity, assuming the loop open at the nonlinearity input. 
If the percentage distortion in this signal is less than 5% then the DF prediction should be valid. It has also been shown 
theoretically that the estimate for the limit cycle frequency, w, is more accurate than for the amplitude, a. Further, as is 


often not clearly pointed out, a is the estimate for the fundamental component in the limit cycle not the peak amplitude. 


When the DF method does not predict a limit cycle it will be correct ifthe C(a) and G(jw) lociare not near to intersecting. 
If this is not the case, however, one can assume a possible limit cycle at the nonlinearity input with amplitude a and 
frequency w corresponding to the near intersection values of C(a) = G(jw). Using the distortion criteria given above a 
very low value must exist for the prediction to be valid when the loci are very close. More details on these concepts are 


presented in section 5.5. 


4.5 Some Examples of Limit Cycle Evaluation 


Some examples are given in the following sections to illustrate limit cycle determination using the DF method for a 
feedback loop containing a single nonlinear element. The first example taken is the Van der Pol equation discussed in 
section 2.2, although a first harmonic balance approach rather than the DF is used as the nonlinear term involves two 
‘inputs, x and x. It is possible, however, to use DF theory for a nonlinear element defined by n(x,x) and the interested 
reader is referred to the first source in the bibliography. The second and third examples consider limit cycles in relay 
systems. It is shown in chapter 6 how limit cycles in relay systems can be calculated accurately so that comparisons of 
the results allows some comments to be made on DF accuracy. The final example is introduced to show the value of the 


IDF in determining limit cycle stability. 


4.5.1. The Van der Pol Equation 


If a harmonic balance approach is used for equation (2.4) one may assume x; = asinwt and x2 = x1 = awcoswt and 


the equation can be written 
X2 +x = pwd = x1)x 


Substituting for x; and x: and neglecting frequencies higher than w in the right hand side gives 


3 
aqw 


—aw sinwt + asinwt = “(aw cos wt — cos wf). 
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For this to be true requires w = 1 for the left hand side and aw = a’ w/a, giving a = 2, for the right hand side. Thus, 
a limit cycle solution, which is independent of y, with frequency 1 rad/s and amplitudes of 2 for x; and x2, is predicted. 


The result is therefore only reasonably accurate for small values of 4, where the phase plane plot is almost elliptical. 


4.5.2 Feedback Loop Containing a Relay with Dead Zone 


For this example the feedback loop of Figure 4.1 is considered with n(x) a relay with dead zone and G(s) = 2/s(s + ty 
. The DF for this relay is obtained from equation (3.20) with A = 0 and is given by 


N(a) = 4h(a’ — &)'" Ja’ x for a > 6, 


which is real because the nonlinearity is single valued. A graph of N(a) against a is given in Figure 4.2. 


n(x) N(a) 


Figure 4.2 Relay with dead zone and its DF 


It shows that N(a) starts at zero, when a = 0, increases to a maximum, with a value of 2h/zd at a = 0 oe. , then 
decreases toward zero for larger inputs. The C(a) locus, shown in Figure 4.3, lies on the negative real axis starting at — oo 
and returning there after reaching a maximum value of — 20/2h. The given transfer function G(jw) crosses the negative 
real axis, as shown in Figure 4.3, at the point (-1, 0) at a frequency of tan"'w = 45° that is w = 1 rad/sec, and therefore 


cuts the C(a) locus twice at the point P. The two possible limit cycle amplitudes at this frequency can be found by solving 


2 
aia i 
2 GG) lo=1 = 1 


Ah(a’ — 6°) 


which assuming 0 = | and h = z gives a = 1.04 and 3.86. Using either the Loeb or the IDF criterion shows the smaller 
amplitude limit cycle to be unstable and the larger one to be stable. If a condition similar to the lower amplitude limit 


cycle is excited in the system, an oscillation will build up and stabilize at the higher amplitude limit cycle. 
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40} Nyquist 


Figure 4.3 Solutions for the two limit cycles at P one unstable the other stable 


The exact frequencies of the limit cycles, obtained using the approach given in Chapter 6, for the smaller and larger 
amplitudes are 0.709 and 0.989 respectively. Although the transfer function is a good low pass filter, the frequency of the 
smaller amplitude (unstable) limit cycle is not predicted accurately because the output from the relay is a waveform with 
narrow pulses, which is highly distorted. If the transfer function of G(s) is K/s(s + 1)’, then no limit cycle will exist in 
the feedback loop, and it will be stable according to the DF method if 


i < Ad (4.7) 
wo(l+o)|,_, 2A 
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that is, K < zd/h. Again, if 6 = 1 and h = z, then the system is stable if K < 1, which may be compared with the exact 
result for stability of K < 0.96. 


4.5.3 Feedback Loop with On Off Relay 


Consider the feedback loop of Figure 3.1 with n(x) an ideal relay, G-(s) = 1, and the plant with a transfer function 
Gi(s) = 10/(s + 1)’. Equation (3.20) with 6 = A = 0 gives N(a) = 4h/az. Thus the C(a) locus, — 1/N(a) =— az/4h 
and lies on the negative real axis starting from the origin with a = 0. The Nyquist locus G(jw) intersects it where it has 


a phase of -180°. The values of a and w at the intersection can be calculated from 


—arl4h = = 7 which can be written 
(1+ jw) 
10 
arg——~— = 180° (4.8) 
"+ jo)? 
and 


an _ 10 
Ah (1 of no ae . 


(4.9) 


The solution for wo from equation (4.8) is tan ' wo = 60°, giving wo = V3 = 1.732. Because the DF solution is 
approximate, the actual measured frequency of oscillation will differ from this value by an amount which will be smaller 
the closer the oscillation is to a sinusoid. The exact frequency of oscillation is 1.708 rads/s so that wo is in error by 1.41%. 
A symmetrical square wave has a harmonic content which is 1/n'* that of the fundamental for n odd. It is thus easy to 
calculate the harmonic components at the output of G(s) and the third and fifth are 5.40% and 1.21% of the fundamental. 


Solving equation (4.9) gives the amplitude of the assumed sinusoidal limit cycle a as Sh/z. 


If the relay has an hysteresis of A, then putting 0 = 0 in equation (3.20) gives 


2 2, 1/2 
N(a) = ante 7 x) gan from which 
an an 
=1 —Tl, 2 Qed: 2 is 
C(a) = —— = — —A + jA. 
OPA a, [ia — A°)'? + jA] 


Thus on the Nyquist plot, C(a) is a line parallel to the real axis at a distance 7A/4h below it, as shown in Figure 4.4 for 
A = landh = z/4 giving C(a) =— a == j. If the same transfer function as above is used for the plant, then 
the limit cycle solution is given by —(a’ — 1)'” —j = im z where w = 1.266, which compares with an exact 


~ (1+ jo) 
solution value of 1.254, and a = 1.91. 
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Nyquist Diagram 


Imaginary Axis 


Real Axis 


Figure 4.4 Nyquist plot for 10/(1+s)? and DF for relay with hysteresis 


For the ideal relay, equation (3.30) gives Niy(a) = 2h/am for the IDF so that the roots of the characteristic equation 
1 + Niy(a)G(s) = 0 show the limit cycle is stable. This agrees with the perturbation approach which also shows that the 


limit cycle is stable when the relay has hysteresis. 


4.5.4 Use of the IDF 


Consider an unstable plant with the transfer function 
G(s) = KI(s — b)(s” + 2s + 10) (4.10) 


in a feedback loop with an ideal relay with output +1. The characteristic equation for the loop with the relay gain taken 


as its DF gives 
s +s (2—b) +s(10 — 2b) +k - 10b = 0 (4.11) 
where KN(a) is replaced by k. Ifa third order characteristic equation has a limit cycle at a frequency , then its form will be 
(s +@)(s +c) =0 (4.12) 
which gives 


(2 — b)(10 — 2b) = k — 10b so that k = 20 — 4b + 2b’. The characteristic equation is stable for b < 2 and the frequency 


of the limit cycle is (10 — 2b)'*. However, when the limit cycle is established the gain through the relay to any unrelated 


signal is the IDF, which for the ideal relay is easily shown from equation (3.30) to be N(a)/2, that is it is halved and the 


new characteristic equation is 


si+s’(2—b) +s5(10—b) +10—-2b4+b’-10b = 0 (4.13) 
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which is stable for b < 0.90. Thus, the odd symmetrical limit cycle cannot exist for b > 0.9. 


Simulation results showed an odd symmetrical limit cycle to around b = 0.4 and an asymmetrical limit cycle for values 
of b > 0.4 up to around b = 0.9, although the limit cycle was dependent on the initial conditions. For b > 0.9 the relay 
locked at either +1 so that the output increased indefinitely. A good explanation for the phenomenon can be seen if the 


root locus for the transfer function G(s) is plotted, as done in Figure 4.5 for K = 1 and b= 1. 


Root Locus 
8 T T T T T T T T T 


System: g | 
6L Gain 17.5 | 
Pole: -0.0268 + 2.621 | 

Damping: 0.0095 
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Frequency (rad/sec): 2.82 | 


System: g 
Gain: 10.2 
“2, Pole: -0.0213 qT 
Damping: 1 = 
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Frequency (rad/sec): 0.0213 | 


Imaginary Axis 
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Figure 4.5 Root locus plot for transfer function of equation (4.7) 
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It can be seen from the gains marked that the complex poles cross the imaginary axis for a gain of 18 and for that value 
of gain the real pole is at -1. However for a gain of 9 (half of 18) the real pole will be in the rhs of the s plane so that the 
system is unstable. If the plot is done for b = 0.9 then the gain when the complex roots cross the axis will be twice that 
for the real root crossing the axis. The difficulty of achieving an odd symmetrical limit cycle for b slightly less than 0.9 is 


presumably due to errors in the DF and the small range of attraction of the limit cycle, which is a closed curve in three space. 


4.6 More than one Nonlinear Element 


When the closed loop system has more than one nonlinear element it is possible to write down nonlinear algebraic 
equations based on the DF method to give any possible limit cycle solutions. The difficulty with this approach, however, is 
that satisfactory accuracy can only be obtained if, when a limit cycle occurs, it is reasonable to assume the inputs at all the 
nonlinear elements are nearly sinusoidal. In the classical multivariable version of Figure 3.1 where G(s) is a transfer matrix 
and n(x) a matrix with nonlinear elements on the diagonal this may be reasonable since typical transfer function elements 
of G(s) will be low pass, which filter the higher harmonic content. However, one has to assume that the limit cycle is at 
a single frequency, which is normally the case in a multivariable control system feedback loop. A simple counterexample 
is easily obtained by lightly coupling two single loops which have limit cycles at two different frequencies. Dependent on 
various factors the loops may continue to maintain both frequencies or lock to a limit cycle with two related, say 1 to 3, 


dominant frequencies. 


Ina single loop configuration the assumption of sinusoidal signals into all the nonlinear elements is usually not justifiable. 
If, for example, another nonlinear element is included in Figure 3.1 prior to G(s) then it may be reasonable to assume 
this has a sinusoidal input but since most controller transfer functions G(s) are not low pass filters the signal into the 
nonlinearity before the plant may be quite highly distorted. Thus, if a DF analysis is done the results may be incorrect. 
A well known example of this is shown in the Simulink diagram of Figure 4.6 where both nonlinear elements are ideal 
saturations with limit levels +1 and the linear elements are a phase lead with a transfer function (1 + s)/(1 + 0.15) 
and a plant with a transfer function of K/s(1 + s)” with K = 9. The linear system and the nonlinear system without the 
saturation characteristic before the phase lead are both stable for K < 11. However the system of Figure 4.6 possesses a 


limit cycle for large initial conditions on the integrator. 


Gain Saturation2 Saturation1 Integrator 


Transfer Fen Transfer Fen 


Figure 4.6 Simulink diagram of system with two ideal saturations 
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Responses for system of Fig 4.6 


Figure 4.7 Initial condition responses of system of Figure 4.6 


Figure 4.7 shows the responses of the system of Figure 4.6 from initial conditions on the output integrator of 4 and 1 
respectively. A limit cycle results from the former whereas for the second case the response is stable although very oscillatory. 
The third response shown on the figure labeled ‘no sat’ is again for an initial condition of 4 but without the saturation 
before the lead network and is again very oscillatory. The limit cycle results because although it is almost sinusoidal at 
the integrator output it is highly distorted at the input to the second saturation, with the result that the fundamental 
component under goes a phase shift when passing through the second saturation. Figure 4.8 shows the input and output 


waveforms for the limit cycle at the second saturation. 


Input and output waveforms of second satuaration 
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Figure 4.8 Waveforms of limit cycle at the input and output of the second saturation 


If with more than one nonlinearity in the loop it is reasonable to assume any limit cycle in the system would be approximately 
sinusoidal at the plant input then, as mentioned in section 1.3, if the plant consists of linear and nonlinear elements a 
frequency response model N(a,w) could be used to model it. This could be done either from mathematical modeling or 
experimental data, where if the distortion is small the output is taken as the signal at the same frequency as the sinusoidal 


input. This frequency dependent DF representation may be used in analysis and design, a topic discussed in section 9.5. 


4.7 Applications of SBDF to find Limit Cycles 


In this section some simple examples are considered to show the use of the SBDF in evaluating limit cycles. The first 
example considers a system with an asymmetrical nonlinearity in an autonomous feedback loop and the second an odd 


symmetrical nonlinearity in a loop with an input which causes a biased operation on the nonlinearity. 
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The procedure is to set up two equations, one for the first harmonic and the other for the bias, or dc, signal. The former 


is always the same, namely 
Na(a,y)GYo) =— 1 (4.14) 


but the latter depends upon where any bias signal enters the loop and the value of G(0). Since for a single valued nonlinearity 


Na(a, 7) is real, the limit cycle frequency is given by 

arg Gjw) =— 180° (4.15) 
which if this frequency is w., gives 

Na(a,y) = 1/1 GGa-) | (4.16) 


The frequency of the limit cycle given by the SBDF solution is thus, unlike the amplitude, unaffected by the bias signal, 


which will not normally be true in practice as the harmonic distortion will be affected by y 


4.7.1. Asymmetrical nonlinearity 


The feedback system of Figure 3.1 is considered with no input G(s) = 2.3/s(s + 1)” and the nonlinearity shown in 
Figure 4.9. It is defined by x + 0.5x” + 0.25x° for |x| < 1, 1 for x > 1 and -1 for x < -1. Figure 4.10 shows the resulting 
limit cycle at the output of G(s) from an output initial condition of 0.8 which since it lies in the range |x| < 1, it may be 


evaluated using its SBDF for this range. 


Figure 4.9 Nonlinear characteristic for Example of section 4.7.2 
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Plot of limit cycle 


Figure 4.10 Limit cycle in Example of section 4.7.1 
For the polynomial nonlinearity n(x) = x” the components of the SBDF from equations (3.26) and (3.27) are 


Nay) = (1) fo + y)"p(xydx 


Nila = ia) f (x + 7)" (ala) p(x) dx 


respectively. After expanding (x + 7)" it can be shown that 
N, (4,7) = 2 nCrte oo 
No(a,7) = (21a") ¥ aCrtters 


where the binomial coefficient 


_ n! 
~ (n—k) tk! 


nC 
Thus for the nonlinearity x + 0.5x° — 0.25x° the values are 
N,(a,7) = 1 + 0.57 + (a’/47) — 3a" /8) — (7°14) 


N.(a,Y) = 1+ y — (0.75a?/4) — 0.757 


(4.16) 


(4.17) 


(4.18) 


(4.19) 


(4.20) 


The linear transfer function has a phase shift of -180° when w = 1 with a corresponding gain of 1.15, so that using equation 


(4.16) gives 


Nila, y) = W1AS = 0.87 


and since in the steady state the average value of the input to the integrator must be zero the bias equation is 


YN, (a,7) = 0 


(4.21) 


(4.22) 
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Equations (4.19) to (4.22) can be solved to obtain a and y. A quick estimate can be found by neglecting terms involving 
powers of y to give y = -0.0836 and a = 0.497. The exact solutions are y = -0.0763 and a = 0.513, which can be seen to 


be in good agreement with the simulation result of Figure 4.10. 


4.7.2 Effect of Bias with Odd Symmetrical Nonlinearity 


The feedback system in the Simulink diagram of Figure 4.11 is considered. The odd symmetrical nonlinear characteristic 
has a slope of 1 for inputs between 0 and 1, 2 between 1 and 2, 1 between 2 and 4 and then zero, and is shown in Figure 
4.12. The linear dynamics are 1.1/s(s + 1)’ and a constant value of 1.5 enters the loop at the nonlinearity output. The 
system is stable when the constant value is zero but its effect is to bias the nonlinearity output so that the average value of 
any limit cycle at the nonlinearity output has a value of 1.5. This is due to the integrator in the loop, the input of which 


must have an average value of zero in the steady state. 


NLfig411 


Constant 


Figure 4.11 Simulink diagram of Example in section 4.7.2 
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Pa: ee es es — ee a 


a 


Figure 4.12 Nonlinearity in the Simulink diagram of Figure 4.11 

The resulting limit cycle waveforms at the input and output of the nonlinearity, when an initial value of 2 is applied at 
the system output, are shown in Figure 4.13. To calculate the simulation result using SBDF analysis is not easy because 
of the linear segmented nonlinear characteristic. It can be simplified using the known simulation solution from which it 
can be seen that the amplitude of the limit cycle at the nonlinearity input is within the range -1 to 2 so that the SBDF just 
needs to be worked out for this range, where the characteristic is a gain of 2 in parallel with minus ideal saturation. The 
bias input is such that it lies on the saturated level of the saturation characteristic, thus y > 1. The bias and fundamental 
outputs can be calculated to be 
bias = 1 — (alx)(1—[(y — 1) /a]’)'” and fund = (2a/z)[(x/4) — (g/2) — (sin 29/4) | 
where sing = [(y —1)/ al. Thus for the nonlinearity the SBDF terms are 

Ny(a,7) = 2 - (ly) + alay) (1 = [7 - Dial)” (4.23) 
and 

N.(a,7) = 2 — (2/x)[(z/4) — (g/2) — (sin29/4)| (4.24) 
G(jw) has a phase shift of 180° for w = 1, with a corresponding gain of 0.55, so that equation (4.16) gives 
Na(a,yv) = 1/0.55 
and the bias equation is 


YN, = (a,7) = 1.5 


The solution to equations (4.23) to (4.26) gives y = 1.20 and a = 0.379 which compares reasonable well with the simulation 


result shown in Figure 4.13 for the nonlinearity input waveform. 
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Limit cycle at the input and output of the nonlinearity. 


Figure 4.13 Simulation result for the example of section 4.7.2 


4.8 Conclusions 


This section has shown how limit cycles can be estimated in nonlinear systems using describing function methods. Systems 
where bias signals may exist in the loop have been evaluated using the SBDF. It has further been shown how the stability 
of any limit cycles can be assessed. Further, it has been demonstrated that when single valued nonlinear elements exist 
in a feedback system failure of the DF method to predict any limit cycles is due to their waveform not being ‘sufficiently 
sinusoidal at a nonlinearity input so that the fundamental signal undergoes a phase shift in passing through the nonlinearity 
which is not predicted by the DE Many textbooks cover some basic uses of the DF usually for evaluating odd symmetrical 
limit cycles and the stability of limit cycles is normally commented on in terms of the Loeb criterion. Some of the more 


detailed topics in this chapter can be found in the books given in the bibliography. 
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5 The SSDF and Harmonically Forced 
Systems 


5.1. Introduction 


Describing functions have so far been considered for a single sinusoid and a sinusoid plus bias. To examine various 
phenomena which may occur in nonlinear systems it is however required to evaluate describing functions for other inputs, 
in particular two sinusoidal inputs. The evaluation of the SSDF for any nonlinearity is considered after introduction of 
the topic by use of a cubic nonlinearity. A general theory is presented but the difficult algebra for specific nonlinearities 
makes evaluation difficult for many nonlinearities, specifically linear segmented ones, so that numerical approaches have 
usually to be adopted to obtain solutions. The problem is shown to be even more complicated if the sinusoidal signals are 
related in frequency, a situation which often needs to be considered in practice. The concept of modified nonlinearities, 
as a means of understanding the evaluation of the SSDF, is then covered and examples given for specific nonlinearities. 
Sometimes it is necessary to consider a perturbation which is related in frequency to an existing sinusoid so the evaluation 
of this IDF is considered next. This is followed by a practical example of the use of the SSDF to evaluate a limit cycle 


using a two harmonic balance. 
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The closed loop frequency response of a nonlinear system is then considered. The most usual situation is for the system to 
have an output with a dominant fundamental component of the same frequency as the input. This situation is examined 
using the DF and it is shown that over some ranges of amplitude and frequency of the input three analytical solutions 
rather than one to the nonlinear equations may be obtained. In practice this results in the so called ‘jump phenomena 
or ‘jump resonance’ in the frequency response. Finally it is shown how its existence can be investigated using the related 
frequency IDF. Various other possible behaviours that may occur, such as the excitation of oscillations and chaotic motion, 
are discussed based on an understanding of the SSDF and some simulation examples given. Brief mention is also made 


of the possible behaviours when a harmonic input is applied to an autonomous system which has a limit cycle. 


Analysis of a specific system using the methods presented in this chapter is not easy unless the techniques are implemented 
in appropriate numerical programs, however, presentation of the theoretical background should allow the reader to 
understand the mechanisms that cause the aforementioned unique behaviours of nonlinear feedback systems. 
5.2 The Cubic Nonlinearity with Two Sinusoidal Inputs 
Consider the two sinusoidal inputs x = acos@ and y = bcos(& +g), where 0. = Wat and 4, = wot into a cubic 
nonlinearity, then the output, f(t), is given by 

f() = la cos 0. + bcos(@ + )} (5.1) 
Expanding the expression yields 


f(t) = a cos’ 0, + 3a’ b cos’ A, cos (8, + ~) + 3ab’ cos O, cos’ (9 + v) + b’ cos’ (8, + 9) 


which on expanding the powers of cos @ may be written 


f(t) = (Bal4)(a° + 2b°) cos A, + (3b/4)(b° + 2a”) cos(@, + ) + (3a”b/2) cos 26. cos (A + @) 

+ (3ab’/2) cos 6, cos 2(0, + @) + (a°/4) cos 30, + (b*/4) cos 3(G, + 9) 

and on using 2 cosA cos B = cos(A + B) + cos(A — B) gives 
f(t) = (Bal4)(a° + 2b”) cos 0, + (3b/4)(b° + 2a”) cos (, + p) + (3a"b/4) 
[cos (26, + 6, + v) + cos(26, — & — ¢)| + (3ab? /4)[ cos (Oa +26, + 20) + cos(O, + 26, + 29)| (5.2) 
+(a°/4) cos 30, + (b°/4) cos 3( + 9) 


If the frequencies w, and w, are unrelated then the first two terms produce the output at the same frequencies as the 
inputs and the SSDF has the components 


N,(a,b) = (3/4)(a’ + 2b”) (5.3) 


N,(a,b) = (3/4)(b’ + 2a’) (5.4) 
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However, if there is an integer relationship between the two frequencies then it is possible that some of the other terms 
in f(t) can produce outputs at the same frequency as the input. For example, if w, = 3. then the second component 
of the third term is (3a°b/4)cos(— 6, — @), which has the same frequency as the sinusoidal input x. Further this 
output component is no longer in phase with the input at frequency w, and the total output at this frequency is 
(3a/4) (a + 2b) cos 0. + (3a’ b/4) cos (Ox + ¢) so that the DF for the input x, which will be denoted N. (a | b) to indicate 


that the signals a and b are related, is 
Na(a|b) = (3/4)(a" + 2b* + abe”) (5.5) 


Thus unlike the unrelated case the above DF value depends not only on the amplitudes of the two sinusoids but also their 
relative phase. It is the phase shift produced through a nonlinearity to an harmonic input related to another harmonic 
input that results in many of the interesting phenomena in nonlinear feedback loops. Further when the nonlinearity is 
not a simple mathematical function, such as the cubic of the above example, it becomes very difficult to get the output 
by the simple trigonometric approach used above. Fortunately, however, it can be shown that for any nonlinearity n(x) 


the output f(f) can be written 


f(t) = > ¥ EjEKQ jx COS JO cos k(O, + ) (5.6) 
j=0 k=0 
where 
a b 
ga f i. ne pT Caley Clb pips ded (5.7) 
—a —b 


and the Neumann factor &, is 1 for n = 0 and 2 for n = 1. 
Comparing equation (5.6) with the result of equation (5.2) for the cubic it can be seen for the cubic that 
aio = (3a/8)(a’ + 2b”), a1 = 3a°b/8,a30 = a°/8 


and the values for @o1, @12, @o3 can be obtained by reversing a and b. These values can be checked to result from equation 
(5.7) with n(x + y) = (x + y)° and it will also be found that Qe = Ofor 7 +k > 3. 


Thus even though the expression of equation (5.6) can in principle be evaluated for any nonlinear characteristic it can 
consist of a large, and possibly infinite number of terms. Also for related input frequencies there may be a large number 


of these terms that produce an output at one of the input frequencies all of which in general will have a different phase. 


5.3. Modified Nonlinearities and the SSDF 


The evaluation of @ from equation (5.7) can be regarded as the two stage process 


a 


n(a,7); = f n(x + 7)Ti(xla)palx)de (5.8) 


—a 


b 


ine f n(a, 7) ;Te(y/b) poy) dy (5.9) 


—b 
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Or, alternatively 


da 


nb = | ny + NROIb»pody (5.10) 


—a 
b 


an = | n(b,yuTi(xla)palx) dx (5.11) 

—b 
The terms n(a, 7); and n(b, vy), may be regarded as generalized modified nonlinearities with the one with subscript zero, 
already having been introduced in section 3.7, using a notation omitting the zero subscript, which will continue to be 
used. Assuming, where unless mentioned to the contrary the two sinusoids are not related, then the SSDF components 
are N,(a,b) = 2@\0/a and N,(a,b) = 2a:/b. An advantage of looking at the derivation in this way is that if one wishes 


to get the SSDF component for x, then one calculates the DF for the original nonlinearity modified by y, namely the DF 


a 


of n(b, y), or alternatively one calculates n(a,y)1 = i; n(x + y)(x/a) pa(x) dx, which is 0.54 times the DF for x in the 


presence of a bias and then averages this with respect to b, as given by equation (5.9) with k = 0. 


The evaluation of modified nonlinearities for linear segmented characteristics is quite difficult and therefore it becomes 
even more difficult to find the SSDF, so that numerical approaches have to be used. However many nonlinear characteristics 
can be approximated over a given range by either a power series or Fourier series expansion, so results are derived below 


for a power law and harmonic nonlinearity, as well as for an ideal relay. 


Download free ebooks at bookboon.com 


71 


An Introduction to Nonlinearity in Control Systems The SSDF and Harmonically Forced Systems 


5.3.1 Power Law Characteristic 


Consider the nonlinearity n(x) = x”. Then one has 
nay = | e+ ¥)" Tala) pas)dx (5.12) 


which on expanding gives 


n(aye= Xf mCax" 7" "Te(xla) pax) dx (5.13) 
n=0 
where the binomial coefficient ,C, = m!/(m — n)!n! Thus the integral can be computed for any chosen value of k, since 


after substituting for T,(x/a), the integral will involve moments of the distribution of p (x). In particular for k = 0, one has 


n(a,7) = x mCnx' Y"" palx)dx = 2 mCnbln (Q) yw" (5.14) 
and for k = 1 
n(a, Yi = » | OP a nea = YC ea) a (5.15) 
n=0 ae n=0 


Using equation (5.14) in equation (5.09) and multiplying by 2/b to get N,(b, 7) from ao: then gives for the SSDF to the 
signal y 


b 
Ns(a,b) = (21b) f¥ mCubtn(ady""B" ypaly)dy = (216?) ¥nCntta( a) ftnnei(B) te) 


—n 


Or, alternatively, equation (5.15) with b replacing a can be used in equation (5.11) and again multiplied by 2/b to get 
N,(b, y) from Qo. to give 


N, (a, b) = (2/b) ri > ed tye 1 (b)x"""Da (x) dx = (2/b°) > mCnn+1 (b) Lm —n(@) (5.17) 
_g 0 n=0 
Since mCn = mCmn-—n these expressions are easily shown to be identical. N, (a,b) can be similarly found or simply obtained 


from N,(a,b) by interchanging a and b. 


5.3.2. Harmonic Nonlinearity 


For the harmonic nonlinearity n(x) = sinmx 


a 


n(a,y)x = fsinmox + y)T(x/a) pa(x) dx 


-a 


which on expanding sinm(x + 7) and substituting x = acos@ gives 


n(a,y)p = ain i cos my sin(ma cos 0) cos k@d@ + i sin my cos (ma cos @) cos co]. 


0 0 
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Substituting the expressions 
cos(zcos @) = Jo(z) + 22 1)‘ Joi(z) cos 210 and 
sin(zcos 0) =— 2>(- 1) fos (z)cos(2i — 1)6 and integrating gives 
n(a, yx = (— 1)*4(ma) sin{ my — (kx/2)} 
In particular for k = 0 the modified nonlinearity is 
n(a,7v) = Jo(ma)sinmy 
which differs only in amplitude by a factor Jo(ma) from the original nonlinearity. Also 
n(a,y); = Ji(ma)cosmy 
It is easily shown, again using the above expansions, that 
N, (a,b) = (2/b)Jo(ma) J, (mb) 


5.3.3 Ideal Relay 


(5.18) 


(5.19) 


(5.20) 


(5.21) 


The two previous characteristics considered were continuous functions of the input variable, whereas this is not the 


case for many approximations to nonlinearities found in control systems, such as the relay or ideal saturation. These 


create difficulties for obtaining general analytical results for the SSDF because the modified nonlinearity will possess 


discontinuities dependent on the relative input magnitudes of the signal and the bias. This is illustrated for the ideal relay 


with output level h for which 


st 


nay = f —hh(xlaypalxddx +f Hli(xlaypalx) dx 
= < 
which gives for k = 1 
n(a,7)x = (2h/kx)Dy(— yla) for |y |< a 
n(a,yY)n = 0 for | |> a 
and for k = 0 


n(a,yv) = (2h/zz) sin | (y/a) for | Y | <a 


n(a,yv) = hsgn y for|y|>a 


(5.22) 


(5.23) 
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If it is assumed that b < a then to find N,(a,b) only n(a,y) is required for | y | < a and 


b b 
N,(a,b) = (2/b) - n(a,y)(y/b) py (y) dy = (4h/7 b’) { ysin '‘(y/a)(b? — yy? dy 
—b —b 


1/2 


Integrating by parts and using the fact that (b” — y’)'” is zero at both limits gives 


1/2 
) 


BD 2 
N,(a,b) = (4h/7° b’) a ane = (8h/ k a){E(k) — (1 — k)K()} (5.24) 
a 
b 


-y) 


where k = b/a < 1 and E(k) and K(k) are elliptic integrals. Similarly, to find N (a,b) only one integral is required for 


b <aifn(a,y), is used and the solution is 


b 2 
N,(a,b) = (4hlaz’) f bo = (8h/az’) E(k) (5.25) 
—b 


Thus even for this very simple break point type nonlinearity the SSDF expressions are rather complicated. 


5.4 The IDF for Related Signals 


The IDF, the gain to an unrelated signal in the presence of a known other one, was considered in section 3.7 for 
the known signal being a sinusoid. Also the perturbation equation of (4.5) to assess stability was derived assuming 
the perturbation was unrelated to the existing signal. Here the IDF is investigated when the small signal is assumed 
to be related to the sinusoid. Consider a nonlinearity with inputs x(t) = acos@, and a small signal d(t). The 
nonlinearity output, making use of a Taylor series expansion, as used to derive equation (4.5), may be written as 
n{ x(t) + d(t)} x= n(x(t)) + d(t)n'[x(o)] where n' (x) = dn(x)/dx. 


The incremental output, dy(t) due to d(t) is thus d(t)n'[x(t)] the output of the nonlinearity n'(x) with input x(t). This 


output may be written for any nonlinearity in the Fourier series. 
n'{(x(t)} = Zc cos (kGu + Px) 
where x(t) has a frequency of @a = Out. If 5(t) is the small sinusoid 0 cos(@, + 6), with w, = Ot, then 
dy(t) = >> cx. cos(k@, + P,) cos(@ + 6) which can be written 
k=0 


dy(t) = 6 3 (cx/2){ cos (kOa + Gc + O + b) + cos(kOa + Gx — 0 — b)} 
k=0 


The only output at w, if @, and Wa, are not related is dco cos(O, + ) but if 61/6, = @a/@, = m/n then the output at 


frequency Wp is 
d{co cos (4, + p) + (cx. /2) cos (6, + Pe - op) | = (5.26) 


which can be written 
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dco cos (O + &) + (Ack /2){ cos (gi — 26) cos (A + b) — sin(Gx — 2g) sin(O + b)}lx=21m } 


Thus, dividing by 6, gives the IDF for situations where the frequencies are related such that if an integer k = 2n/m exists, 


then it is given by 
Nia (a) = Co + (cx/2) es | k=2n/m (5.27) 


where the subscript 6 has been used to distinguish it from the unrelated IDF N,,(). Since in general the coefficients c, 
normally decrease with increasing k, the most important ones are for low values of k. Those for k = 1 and 2, in particular, 
will be of some interest later. From equation (5.27) it can be seen that N,(a) is complex and its locus is a circle centred at 
c, with radius c,/2. If one wishes to plot Cis(a) =— 1/Nis(a) on a Nyquist plot then it can be shown that Cis(a) is also a 
circle with centre at — col {co + (cx/2)°} and radius (cx/2){ c6 + (cx/2)?}. 


Also for a single valued nonlinearity with k = 2, which is the situation for a perturbation of the same frequency, @, = 0 and 
2 =2 - n' (x)[(2x" /a’) — 1]p(x)dx = 2[.Ni,(a) -— N(@)] 
so that 


Ns(a) = N(a) + (a/2)N'(a)(1 + &*) (5.28) 


This is the equation of a circle with centre at (Nj,(a),0) and radius (a/2)N' (a). 
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5.5 More Accurate Determination of Limit Cycles 


The classical approach is to balance more harmonics not just the fundamental. Although straightforward in principle it is 


not easy analytically due to the difficulty of computing the output from a nonlinear element subject to two or more related 


harmonic inputs. To illustrate the procedure a feedback loop with the simple polynomial nonlinearity n(x) = x — (x°/6) 


and a loop transfer function G(s) = 2(1 — s)/s(1 +s) is considered. 
The nonlinearity input is taken as 
x = acoswt + bcos(3at + ¢@) 
Computing the output from the nonlinearity at frequencies w and 3w, gives the terms 
la — (a/8) (a + 2b° + abcos 9| cos wt + (a’b/8) sin @ sin wt 
and 
[b — (b/8)(b? + 2a’) — (a*/24) cos py |cos (3t + g) — (a*/24)sing sin (3at + 9) 
Thus the magnitude and phase components for the SSDF are 
M.(a|b) = {[1 — (1/8) (a? + 2b’ + abcos¢) | +[(1/8)ab sin pf} 
M,(a|b) = (1/b){[b — (6/8) (b? + 2a’) + (a°/24) cos gf + [(a?/24)sin ef} 
dala | b) = tan”! |—ab sin g/(8 — a’ — 2b’ — abcos 9) | 


hy (a | b) = tan”! [a’ sin g/ (24b — 3b*° — 6a’b — a’ cos y)| 


(5.29) 


(5.30) 


(5.31) 


(5.32) 


(5.33) 


(5.34) 


(5.35) 


If the gain and phase of the transfer function are denoted by g and 8 then the harmonic balance loop equations for the 


two frequencies are 
M.(a|b)g(w) = 1 
M,(a|b)g(3o) = 1 
ba(a|b) + O(w) =— 180° 


bv(a|b) + O(3w) = (2j + 1) 180° 


(5.36) 


(5.37) 


(5.38) 


(5.39) 


Download free ebooks at bookboon.com 


76 


An Introduction to Nonlinearity in Control Systems The SSDF and Harmonically Forced Systems 


An estimate for the solution can be obtained by using the DF solution of a = 2 and w = 1, together with equations 
(5.37) and (5.39) with (5.33) and (5.35) substituted. This gives b = 0.222 and @ = -53° or 127°. Solution of the four 
nonlinear equations (5.36) to (5.39) with the values from equations (5.32) to (5.35) substituted gives two solutions of 
a = 2.19, b = 0.298, @ = 138°, w = 0.882 and a = 1.88, b = 0.166, p = -52°, w = 1.06. The former solution is that of the 
stable limit cycle shown in the simulation result of Figure 5.1. Also shown are the simulation results for the output of the 
nonlinearity and the calculated result for the stable limit cycle marked theoretical. Note the phase origin is arbitrary so 


the waveform has been plotted just to show the good waveform agreement. 


Simulation result for two harmonic balance Fig 5.1 


Figure 5.1 Comparison of theoretical and simulation limit cycle results. 


It is relatively easy to write numerical routines to compute more precisely any limit cycles which may exist in a feedback 
loop. References 5.1 to 5.3 give details of some papers which discuss possible approaches using programs which make use 
of the fast Fourier transform and interactive graphics. Because a large number of nonlinear algebraic equations need to be 
solve an initial guess for the limit cycle is required. This can be obtained using the DE, if a predicted solution exists, or if 
not by assuming a solution using say an amplitude and frequency near to an intersection of the Nyquist plot and the C(a) 
locus. The assumed limit cycle solution, x,(t), is then used as input to the nonlinearity with the loop assumed open and 
the signal fed back to this point x,(t) calculated. This is done by discretising the signal, finding the corresponding output 
from the nonlinearity, evaluating its Fourier series, passing these components through G(s) and recombining them to 
produce the nonlinearity input. Graphs of x,(t) and x,(t) are then plotted for comparative purposes. Further iteration can 
be done to see if the waveforms converge, typically they may become quite similar but because the assumed frequency is 


incorrect appear to have a phase difference. 
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The frequency can then be changed to see if convergence occurs. Alternatively at any point in the process a harmonic 
balance may be computed using a limited number of harmonics in x,(t) by taking it equal to the last value of x,(t). Using 
these techniques examples are given in the papers referenced of deriving limit cycles in a large number of feedback systems 
some with multiple nonlinear and linear elements. Cases are considered when the DF method does and does not predict a 
limit cycle. An interesting feature of the approach is that because it is numerical it can obtain unstable limit cycle solutions 
which cannot be found by simulation. One example which was considered is the feedback loop originally investigated by 


Fitts with a cubic nonlinearity and linear dynamics 
G(s) = 10s°/(s* + 0.04s° + 2.0206s° + 0.0404 + 0.9803). (5.40) 


This transfer function, which is not typical of one found for a control system with the double derivative in the numerator, 
has two pairs of lightly damped poles near to frequencies of 0.9 and 1.1 rads/s. and its Nyquist plot is shown in Figure 
5.2. The feature of this Nyquist plot, which causes the DF method to fail, is that it has related frequencies on different 
sides of the C(a) plot, which for a cubic lies on the negative real axis. The Nyquist plot starts from the origin and loops 
in the positive part of the plane to the real axis which it reaches at a frequency of 1.005 before looping in the negative 


half plane to finish at the origin. 
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Figure 5.2 Nyquist plot for transfer function of equation (5.40) 


In chapter 6 it is shown how limit cycles can be calculated exactly in feedback loops with a relay nonlinear element. Any 
nonlinear characteristic can be approximated by a quantization characteristic, the more quantization levels used the more 
accurate the approximation. Since a quantization characteristic can be obtained by relays in parallel the method for finding 
limit cycles in relay systems can be used for any nonlinearity using a quantized approximation. The problem with the 
method is that the number of nonlinear equations which have to be solved increases with the number of quantization 
levels used. However, by slowly increasing the number of quantization levels good initial estimates can be obtained for 


the nonlinear equations and good results are reported in reference 5.4 using this approach. 


5.6 Closed Loop Frequency Response 


When a sinusoidal input is applied to a feedback loop which is linear then provided it is stable the output will also be a 
sinusoid with amplitude and phase dependent on the input frequency. Further if the input amplitude is doubled then this 
just doubles the output amplitude. Obviously if the system is nonlinear then the frequency response will be amplitude 
dependent but an obvious question is will the output be roughly sinusoidal and at the same frequency as the input? This 
aspect will be examined in the next section, 5.7,by assuming that it will be the case but before doing so it is appropriate 


to comment on some other possibilities suggested by the previously derived DF results. 


Results for the SSDF have shown that the both the gain and phase shift which a sinusoid undergoes in passing through a 
nonlinear element can be affected by another signal. If the other signal is unrelated to an existing one then only the gain to 
the existing one will be affected For a limit cycle to build up in a feedback loop then the loop must be unstable to a small 
signal, clearly a situation which might arise when an external signal causes the gain through a nonlinear element to a small 
signal to change appropriately. Figure 5.3 shows a simple feedback loop with a nonlinear element n(x) = x° for |x| <1, 1 


for |x| > 1 and -1 for |x| < -1, dynamics of 2/s(s + 1)’ and a sinusoidal input of amplitude 1 at a frequency of 4.1 rad/s. 
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Sine Wave Integrator 


Transfer Fen Scope 


NLcubsat 


Figure 5.3 Simple feedback loop with cubsat nonlinearity 


Although the system is stable with no input, as indicated by the DF method, or as can be ascertained by simulation, when 
subject to the above sinusoidal input the gain through the nonlinearity is increased sufficiently to cause a limit cycle to 
occur at an approximate frequency of 1 rad/s, the 180° phase shift frequency of G(s). The output of the system showing 


the dominant limit cycle frequency, which is unrelated to that of the input, is shown in Figure 5.4 


Output from system of Fig 5.3 
r r — 


—— 


Figure 5.4 Output of System of Figure 5.3 


Also the value of Nis(a) for k = 1, which is achieved with 2n = m, can have a phase shift. This indicates the possibility of 
a limit cycle occurring which is at one half of the frequency of the applied input, namely a second subharmonic. Clearly 
also other subharmonics might exist if they can somehow be excited since a phase change is possible for finite amplitudes 


of other related frequencies, such as the ratio 3:1. 


A further behaviour that is possible is for the output to have a chaotic motion, that is a waveform which is deterministic, 
in the sense that if the experiment is repeated the same waveform will result, but it cannot be described analytically. This is 
demonstrated for the system of Figure 5.5 with an unusual nonlinearity n(x) = x — (x°/6), which changes from a positive 
to a negative slope. The loop dynamics are 2/(1 + s)° and the input frequency of the sinusoid is 1 rad/s. Figures 5.6 to 
5.8 show the output response as the amplitude of the input sinusoid is increased. In Figure 5.6 for an input amplitude of 


1.1 the output, unlike for smaller input amplitudes, is noticeably distorted from a sinusoid. 
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Figure 5.5 System which demonstrates chaotic motion. 
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Figure 5.6 output for input amplitude of 1.1. 
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With an increase in amplitude to 1.3, shown in Figure 5.7, the output has reached a bias level and shows considerable 
third harmonic distortion. Finally in Figure 5.8 for an input amplitude of 1.4 the output shows a chaotic motion and for 


an input amplitude of around 1.5 the response goes unbounded. 


Forced response 


Amplitude 1.3 


Figure 5.7 Output for input amplitude of 1.3 


Frequency Response 


Amplitude 1.4 


Figure 5.8 Chaotic motion for an input amplitude of 1.4. 


Obviously another possibility is that one may wish to apply a sinusoidal input to a system that already has a limit cycle, a 
situation which once more means two sinusoidal inputs to a nonlinearity need to be considered. Again based on the SSDF 
three possibilities can be expected namely (i) the limit cycle and forcing signal will coexist, (ii) the forcing signal may 


‘quench’ the limit cycle or (iii) the forcing signal and limit cycle may ‘lock to related frequencies, known as synchronization. 
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Clearly to investigate theoretically any of these varied behaviours which may exist in a nonlinear system for a sinusoidal 
input one has to apply the appropriate analysis to see if such a behaviour might exist. Secondly one has then to show 
that the behaviour is stable. Affirmative solutions to these two questions then means, assuming the approximate analysis 
is sufficiently accurate, that the behaviour will probably exist. However the analysis is based on steady state assumptions 
and whether a specific steady state is reached is a dynamic problem which is difficult to solve. Like singular points more 
than one behaviour may be reached dependent on the initial conditions. Luckily most of the above behaviours are rare 
in a feedback control system since to exist they require features such as a lightly damped feedback loop, or ‘unusual’ 
nonlinearities, for example possessing both positive and negative slopes. In the next section the ‘normal frequency response 


is examined such as that found for the system of Figure 5.5 for input amplitudes less than unity. 


5.7. Jump Resonance 


When the closed loop system of Figure 3.1 has a sinusoidal input r(t) = Rsin(wt + f), it is possible to evaluate the closed 
loop frequency response using the DF. If the feedback loop has no limit cycle when r(t) = 0 and, in addition, the sinusoidal 
input r(t) does not induce a limit cycle, then, provided that G.(s)Gi(s) gives good filtering, x(t), the nonlinearity input 
may be taken equal to the sinusoid a sin wt. Thus with the nonlinear element replaced by its DE, equations can be obtained 


for the loop operation. For example balancing the components of frequency w around the loop gives 


Reg-sin(wt + 8 + &) — aM(a)gigcsin(wt + f(a) + & + @) = asinwt (5.41) 
where G, (jw) = gue and G,|(j@) = g1 e””. Then, in principle, equation (5.41) which can be written as two nonlinear 
algebraic equations, by resolving in the directions of wt and cos wt, can be solved for the two unknowns a and f for 


the given sinusoidal input with magnitude, R, and frequency, w. The fundamental output signal can then be found from 
c(t) = aM(a)gisin(wt + d(a) + A) (5.42) 


Various graphical procedures [5.5-5.7] have been proposed for solving the two nonlinear algebraic equations resulting 
from equation (5.41) or for the loop balance results written in a slightly different form. If the system is lightly damped, 
the nonlinear equations may have more than one solution, indicating that the frequency response of the system may 
have a jump resonance. This phenomenon of a nonlinear system has been studied by many authors, both theoretically 
and for possible practical applications. The major advantage of a graphical approach to obtain the solution, rather than 
a computational approach to solving the nonlinear algebraic equations of (5.41), is the understanding it provides of the 
jump resonance phenomenon. One approach is briefly outlined which is applicable for a single valued nonlinearity, that 


is for 6(a) = 0. From equation (5.41), on separating the sinwt and cos wt terms, one obtains 
Rg.cos(B + 6-) — aM(a) gi: g-cos(( + 6-) = a (5.43) 
and 


Re-sin(S + @) = aM(a)gig-sin(a + &) (5.44) 
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Denoting the normalized input magnitude to the nonlinearity a/R by x, the normalized output aM(a)/R by y and 


eliminating f gives 
i= ge[1 —y gi sin’ (A + 6.) |? — ygig-cos(@ + @.) (5.45) 


This is the linear frequency dependent relationship between the nonlinearity input and output which needs to be solved 


together with the DF relationship 
y = xM(xk) (5.46) 


which needs to be plotted on the same x-y plot as equation (5.44) for the given value of R to get possible solutions. In 
cases where jump resonance occurs there are typically three solutions rather than one solution [5.8-5.10]. In this case 
it is possible using the stability criterion of the next section to show that two solutions will be stable and one unstable. 
An interesting phenomenon called island-type jump resonance is described in reference 5.11. Here one stable solution is 
isolated, so that it may only be reached from the application of suitable initial conditions, and therefore may not be seen 


in a standard frequency response test. 
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5.7.1 Jump resonance region- use of the IDF 


Determining whether a specific system can have a jump resonance in its frequency response can be quite laborious using 
the method of the previous section. Based on a phasor diagram for the DF linearised closed loop system it can be argued 
that for a constant frequency input if its magnitude R is increased then a should also increase. Thus Hatanaka [5.12] 


obtained the criterion that for a system to possess jump resonance 


(sr 


facts but all are consistent with the fact that jump resonance will occur ifa solution exists for the IDF characteristic equation 


wconst 


R 
) < 0 and the jump resonance points will be given by (o%) =0. Further results have been found using these 
q@ceonst 


1+ Nis(a)GUiw) = 0 (5.46) 
where the IDF is for a perturbation of the same frequency as the input as given in equation (5.28). 


For a specific solution amplitude a this equation is easily checked but of course this requires the solution for a to be 
available. By plotting Cis(a) =— 1/Nis(a), which have been shown to be circles, for all possible values of a then no 
intersection of G(jw) with these loci ensures that jump resonance will not exist. From equation (5.28) the maximum 


phase shift produced for a given a is @ where 


| (a/2)N' (a) | 


(5.47) 
N(a) + (a/2)N' (a) 


sing = 


If @m is the maximum value of @ for all a then no jump resonance will exist if G(w) does not enter the sector on the 
Nyquist plane bounded by the radial lines 180 + @m. Interestingly for a power law nonlinearity n(x) = c| x |“sgn(x) it 


can be shown that 
a = tan '[(1 — w)/2p'7] (5.48) 
which is independent of a so that @m = @ for all a. In particular for a cubic “4 = 3 and a = 30°. 


5.7.2 Some Examples of Jump Resonance 


As stated above the DF analysis presented in the previous sections is based on steady state conditions. Whether a specific 
steady state is reachable is a dynamic problem. The existence of jump resonance behaviour was reported in the 50’ and 60’s 
from some studies, amongst others, of nonlinear electrical circuits and analogue computer simulations. The experiments 
were usually done using sinusoidal signal generators with the amplitude and frequency controlled manually. Typically 
the amplitude for the sinusoidal signal generator was set at a fixed value and the frequency varied slowly by hand turning 
of the frequency setting dial, first in an increasing direction to a certain value and then decreasing back to the starting 
frequency. Clearly there was no careful control of the dynamics of the procedure but two different levels of output were 
found over the ‘jump range; one for increasing and the other for decreasing the frequency. Digital simulation allows more 
precise control of the dynamics of an experiment so for investigating the jump phenomenon questions arise as to how the 


sweep should be generated and how fast it should be to ‘claim’ the steady state has been reached. 
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I am indebted to H.K.Wong of the University of Warwick for the majority of the results in this section. To start the study 
a simple feedback loop having a transfer function G(s) = 1/s(s + c) in the forward path and a cubic nonlinearity in the 


feedback path is considered. A simulink diagram is shown in Figure 5.9 with c = 0.4. 


1 
s+0.4 


Sine Wave Transfer Fen Integrator 


Scope 


Fen To Workspace 


Figure 5.9 Simulink Block Diagram for Simulation Study 


Denoting the input to the NL by acos wt, which is the system output, and the input to the system by Rcos(wt + £), the 


loop balance equation is 


Re cos(wt + 8 + 6) — aM(a)gcos(wt + @) = acoswt (5.49) 
giving 

Rg cos(B + 0) — aM(a)gcos 6 = a (5.50) 
and 

—Rsin(B + 6) + aM(a)sin@ = 0. (5.51) 


jO(o) 


where G(jw) = g(w)e’” and the dependence of g and @ on w has been omitted. 


Denoting aM(a)/R as y and a/R as x and eliminating f from the above equations gives 
y? + (x7 /g*) — (Qxy cos O/g) = 1 (5.52) 


To show graphically the solutions this x-y relationship for the frequency response of G is plotted for different values of w 
and the x-y relationship for the cubic is then superimposed on it. The results are shown in Figure 5.10 for a range of w 
together with the cubic relationship for R=1, which is simply y = x * M(x) and is y = 0.75x° for the cubic. By plotting 
the relationship in this way the possible solutions are easily seen for a fixed R and if R is changed only the single curve 


of the nonlinear relationship needs to be changed. 
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Implicit elliptical curves as loci of frequencies representing system state of a 2-element closed-loop 
feedback system consisting of of a G(s)=1/s°+0.4s linearity and a cubic nonlinearity with the latter in the feedback path 
This is overlayed with sinusoidal single input describing function of the nonlineary in its I/O space 


45 n the system is excited by a sinusoidal input of fixed amplitude R = 1unit(s 


—— 0.500 rad/s 
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System output amplitude = NL invut amplitude x 


(Jump region) 


45 
NL output amplitude v 


Figure 5.10 Ellipses for G(j) with c = 0.4 and cubic nonlinearity for system of Figure 5.9 
The intersections between each ellipse for a particular frequency with the nonlinearity input-output (x-y) plot mark the 
amplitude solutions for the fundamental in this x-y space. The coordinates of the intersections can then be used to compute 
the magnitude or phase response of the system. Jump phenomena may only occur in the frequency range where there is 
more than one intersection of the nonlinear characteristic with an ellipse for that frequency 
Substituting M(a) = 3a’ /4 for the cubic in equations (5.50) and (5.51) gives 
Rg cos(B + 6) — 0.75a° gcos@ = a (5.53) 


and 


—Rsin(B + 6) + 0.75a’ sind = 0 (5.54) 
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These equations can be solved numerically for a and f for given R and w avoiding the graphical approach. However, 
initial guesses for the solution values have to be given to a nonlinear algebraic equation solver, which is greatly facilitated 
by knowledge from the graphical approach. From software written using the graphical approach [5.13] the theoretical 


frequency response predictions are given in Figure 5.11 for c = 0.2 in the transfer function G and input R = 1. 


Prediction of frequency response for a closed-loop feedback system wth known linear and nonlinear elements: 
; = _ As*+02s, and a cubic nonlinearity with Nonlinearity as feedback _ ; / 
using a combination of sinusoidal singe input describing function and graphical methods for an input amplitude of 1unit(s) 


Frequency (rad/s) 


0 CE 


Phase (deg) 


0 05 1 
Frequency (rad/s) 


Figure 5.11 Solutions from the graphical method for the frequency response with c = 0.2 


A software tool [5.14] was used to obtain sweep frequency response tests in simulink. The tool was modified so that it 
could perform bi-directional sweeps necessary for the jump phenomenon investigation. The input sweep generator changes 
the frequency of the sinusoidal input from @ to @, rads/s progressively so that a =r '@x-1,7 | > 1 for sweep up, and 
from @, to @; so that W.—1 = r@x,r < 1 for sweep down. The frequency increment is only applied after the output satisfies 
a steady state test. This is done by checking over the last n cycles of the output and confirming that the phase change is 
less than a specified tolerance, e. Averaging is also performed for m specified cycles of the output to compute the result. 
To obtain the results n and m were set to 3, € to 10-5 and r to 0.995. When the increment in frequency is performed the 
system state conditions are preserved between frequencies, thus ensuring that conditions are similar to the previously 
mentioned analogue computer approach. The measured output is the fundamental amplitude of the output waveform 


determined by Fourier analysis, which should give best agreement with the DF analysis. 
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The results obtained for the cases of c = 0.4 and 0.2 with R =1 are shown in Figures 5.12 and 5.13, superimposed on the 
theoretical curves. It can be seen that there is good agreement between the theoretical and experimental results particularly 


around the jump frequencies. 


Prediction of frequency response for a clased-loop feedback system with known linear and nonlinear elements: 
A s*+04s, and a cubic nonlinearity with Nonlinearity as feedback 
using a Combination of sinusoidal singe input describing function and graphical methods for an input amplitude of 1unit(s) 


Frequency (rads) 


Figure 5.12 Swept frequency response test for c = 0.4 compared with the theory. 
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Prediction of frequency response for a closed-loop feedback system with known linear and nonlinear elements: 
A s*+0.2s, and a cubic nonlinearity with Nonlinearity as feedback 
using a combination of sinusoidal singe input describing function and graphical methods for an input amplitude of 1unit(s) 


Frequency (rad! s) 


Figure 5.13 Swept frequency response test for c = 0.2 compared with the theory. 


In addition spot tests at various frequencies within the predicted jump range were undertaken. This was done in simulink 
by defining the input as a unit amplitude sine wave with phase zero and simulating G(s) as in Figure 5.9, which allows initial 
conditions to be applied to the integrator. For c = 0.4 two almost sinusoidal output amplitude solutions in good agreement 
with the theory were obtained in the frequency range 1.38 to 1.49 as shown in Figure 5.14 for a frequency of 1.45. The 


smaller amplitude solution was obtained with a zero initial condition and the larger one with an initial condition of -2. 
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“0 5 10 15 20 25 30 35 40 


Figure 5.14 Plots of the two outputs for c = 0.4 and input frequency of 1.45 rads/s. 


For c = 0.2 the larger output amplitude solution could only be found with an initial condition on the integrator for 
frequencies above 1.54 up to 2.1 rads/s. With a zero initial condition on the integrator the output was nearly sinusoidal 
at the input frequency until a frequency just exceeding 1.66 rad/s was reached when the lower amplitude contained a 
significant third harmonic which remained until a frequency of around 1.75 rads/s. Between approximately 1.75 rads/s 
and 1.80 rads/s the third subharmonic disappeared and the output again became almost sinusoidal at the input frequency. 


The third subharmonic component again returned from just above 1.80 rads/s until it finally disappeared at 2.05 rads/s. 


What do you see? 


A skateboarder? 


We see an All Options Junior Trader leaving 
the office after a successful day. 
He has the freedom to make his own decisions. 
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Further investigation for frequencies between between 1.66 rads/s and 2.10 rads/s using different initial conditions on 
the integrator revealed three stable modes according to the initial conditions. These were two approximately sinusoidal 
amplitudes at the input frequency in good agreement with the theory and a third with a significant subharmonic content. 


Figure 5.15 shows these three outputs for an input frequency of 1.9 rads/s. 


Figure 5.15 Plots of the three outputs for c = 0.2, R= 1 and an input frequency of 1.9 rads/s. 


This example clearly illustrates that use of the approximate DF method can determine steady state frequency response 
solutions but cannot predict what can occur in a particular situation, as this is determined by the transient dynamics 
when more than one stable state exists. The possibility of the third solution existing, namely the subharmonic oscillation 


can also be predicted by the DF method as is shown in the following. 
If the input to the feedback loop is again taken as Rcos(wt + f), then the input to the cubic nonlinearity, which is also 
the output, may be taken as acoswt + bcos((wt/3) + @) to investigate the possibility of a third subharmonic. The 
fundamental output signal from the cubic is aM. (a | b)cos(wt + dala | b)) and the loop equation for the fundamental, 
with the dependence of g and 6, on w shown explicitly in the equation, is 

Re(w)cos(wt + 8 + 6(w)) — aMa(a | b) g(w) cos(wt + O(w) + dala | b)) = acoswt (5.55) 
For the third subharmonic the corresponding loop equation is 


—bM,(a | b) g(@/3) cos ((wt/3) + A(w/3) + do(a | b) + v) = bcos((wt/3) 9) (5.56) 


These two equations provide four nonlinear equations which can be solved for the unknowns a, b, § and ¢ to obtain a 
possible subharmonic solution when R and w are given. A quick test for a possible subharmonic can be obtained from 


the phase shift requirement of equation (5.56), namely 


O(w/3) + dv(a| b) =— 180° (5.57) 
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From equation (5.2) it can be shown that 
hy (a | b) =— tan"! [ab sin 39/ (b> + 2a’ + abcos 39) | (5.58) 


Differentiation of the expression in the brackets shows that it isa maximum when cos 3g = — ab/ (b° + 2a’) which gives 
the maximum value for é,(a | b) of +sin ' | ab/ (b° + 2a’)]. Denoting the ratio of b/a by p gives sin | [o/ (0° + 2)], 


which on differentiation gives a maximum when 0 = J? and a maximum value of sin’! (2. /4) = 20.7°. 
8 


Thus a subharmonic can only occur provided the phase lag of G(jw) is greater than ~159°. The lowest frequency for which 
this occurs is approximately 0.52 rads/s, so that subharmonic oscillations might be possible when the input frequency 
exceeds 1.56 rads/s which is in reasonable agreement with the experimental results. The subharmonic only lasts for a 
small frequency range as the loop gain to the subharmonic of M,(a | b) g(w/3) drops below unity. It can be shown that 
the inverse of the describing function for the subharmonic is a set of circles dependent on p within the sector of angle 
+20.7° about the negative real axis of the Nyquist plane. Thus for a subharmonic to exist for the cubic nonlinearity in a 


feedback loop the Nyquist plot of GGjw) must enter this sector. 


The disagreement in Figures 5.12 and 5.13 at low frequencies is due to the poor filtering of the harmonics provided by 
G(jw). If the cubic is moved to the forward loop in front of G(s) the discrepancy at low frequencies is less as the sinusoidal 
loop input signal passes straight to the cubic nonlinearity. This is illustrated by the frequency response shown in Figure 5.16 
for the case of c = 0.4 and R = 1. Here, however, another situation arises in that theoretically there are only two solutions 
above a frequency of around 2.25 rads/s. and thus there is no jump down on the upward frequency sweep as confirmed 
by the sweep frequency test, due to the loop having a negligible feedback signal. In practice, of course the increasing gain 
of the cubic is not possible and a jump down will result. Spot testing shows that if no initial conditions are placed on the 
integrator the solutions of the decreasing frequency sweep are obtained. However above the frequency of 2.25 rads/s the 


high amplitude stable solution can be obtained by putting initial conditions on the integrator. 
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Prediction of frequency response for a closed-loop feedback system with known linear and nonlinear elements: 


A s*+0.4s, and a cubic nonlinearity with a Hammerstein forward path configuration 
using a combination of sinusoidal single input describing function and gaphical methods for an input amplitude of tunit(s) 


Magnitude 


8 ——s 
*® OF prediction 
Sweep up from f=1 
————— Sweep down from f=5 


Frequency (rad/s) 


Figure 5.16 Frequency response for cubic in the forward path with c = 0.4. 
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Results are given in Figure 5.17 for an ideal saturation characteristic in the forward path before the same G(s) with 
c = 0.2 and R = 1. In this case it is seen that the theoretical resonant peak bends in a backward direction so that the smaller 


jump occurs on the upward frequency sweep. 


Figure 5.17 Frequency response results for saturation in the forward path, c = 0.2 and R= 1. 


5.8 Conclusions 


In this chapter the response of a nonlinearity to two sinusoidal inputs was first discussed and although it was shown that 
a mathematical formulation is straightforward obtaining analytical solutions for other than simply defined mathematical 
characteristics is not. Thus numerical computation approaches have to be adopted for most practical situations. The 
concepts were used, however, to show how various properties of feedback systems could be investigated by considering 
two frequencies present in the loop. The analytical work has been supported by simulation results and the chapter was 


concluded by a detailed studied of jump resonance. 
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6 Limit cycles in relay systems 


6.1 Introduction 


In the early days of control engineering, when power electronics was primitive compared with the situation today, relays 
were commonly used as they provided a cheap form of power amplification. They are still used in temperature control, 
where the temperature is controlled to limit cycle around the desired value, and switching techniques using modern power 
electronics are also frequently used. The distinction between switching nonlinearities and other nonlinearities is that their 
input does not control the output continuously but only determines the instant of switching. The input has no control of 
the output between switching instants This feature, that the input has no control of the output between switching instants, 


allows a unique theory to be developed for analysing limit cycles in feedback loops with relay elements. 


The starting point for the analysis is to assume a specific output waveform from the relay and develop equations which must 
be satisfied for the assumed waveform to exist. The most common situation is the case of a relay with no dead zone where 
the output will typically be assumed to be a square waveform with a one to one mark space ratio. More complex modes 
with multiple pulses per cycle can exist in some systems and their prediction is still possible, however more nonlinear 
equations need to be satisfied, thus typically requiring a computational approach. The analysis starts by assuming a periodic 
output waveform from the relay and then the requirement is to obtain the steady state periodic output this produces after 


passage through the transfer function G(s). This is also a problem of interest in circuit theory [6.1]. 


life- 
changing 
careers 


Where your talents benefit you, and millions of others 


At Novo Nordisk, we are working toward a future where everyone's potential can 
be fulfilled. What about yours? Right now, we are looking for outstanding candi- 
dates to join our 7 Graduate Programmes. This is your opportunity to dream big, 
realise your potential and help defeat diabetes worldwide. 


7 International Graduate Programmes 

Graduate Programmes are available within Global Finance, Global Marketing, 
Supply Chain, Business IT, Research & Development, Business Processes and 
Pharma Management Education. Most programmes are based in Denmark and 
Switzerland with rotations abroad. All programmes start 1 September, 2009. 


Read more about our Graduate Programmes at novonordisk.com/graduates. 
Deadline for application is 12 February, 2009. 


www.novonordisk.com/graduates @ 


novo nordisk” 


Download free ebooks at bookboon.com 


97 


An Introduction to Nonlinearity in Control Systems Limit cycles in relay systems 


For application to the relay feedback problem three different methods have been proposed for doing the analysis all of 
which produce identical results. A time domain approach was developed by Hamel [6.2], a frequency domain approach 
by Tsypkin [6.3], and a time domain approach using a state space representation by Chung [6.4, 6.5]. The frequency 
domain approach uses a ‘summed’ frequency locus which makes it easy to compare with the DF method so that this will 
be the method considered below. The formulation is done in terms of A loci which are more general than the Tsypkin 
loci, and this is followed with details on the properties of A loci and their evaluation. Solutions for limit cycles of the 
‘fundamental’ form are discussed. In section 6.5 it is shown that an exact condition can be found giving necessary and 
sufficient conditions for the stability of any predicted limit cycle, which is an extremely useful result. This enables the study 
in section 6.6 of some very interesting problems relating to limit cycles in relay systems which have been obtained from 
software implementations of the theory presented. Topics discussed include limit cycles with multiple pulses per period, 
limit cycles with sliding and concepts for the prediction of chaotic motion. Finally in section 6.7 the harmonically forced 


response of relay systems is considered. 


6.2 The Frequency Domain Approach 


Consider a periodic pulse waveform for which a one period example is illustrated in Figure 6.1, where it is assumed to be 
odd symmetrical and to have eight pulses in a period. If any one of the eight pulses is denoted by y;(t) which is assumed 
to be of width Aj; and to start at time ¢; then it is straightforward to show that it has the Fourier series 


yi(t) = (AAG IT) + (hi) & (/n)[sinnwAt;cos nw (t — ti) + (1 — cosnwAt) sinna (t — t))| (6.1) 
n=1 


If y;(¢) is the input to a linear system with the transfer function G(s), then the Fourier series for the output c;(t), assuming 


lims...G(s) = 0 and G(s) is finite, is given by 
ci(t) = (hAtG(0)/T) + (Al) 

Dy (gn/n){sin nwAt;cos[ nw (t — t) + gn] + ( — cosnwAt) sin[ nw (t — t) + ¢,|} (6.2) 
where the frequency response of G(s) is denoted by 


G(jnw) = gne’™ = Uc(nw) + jVc(nw) (6.3) 


Figure 6.1 One period of an odd symmetrical multipulse waveform. 
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If equation (6.2) is differentiated then the Fourier series for ci(t), provided lim;_..sG(s) = 0, is 
ci(t) = (wh/7) 

p> (2n In){- sinnwlt; sin| nw (t —t)+ n| + (1 — cosnwAt;) cos| nw (t —t) + n|} (6.4) 
When lim;_~sG(s) 4 0, c(t), will have discontinuities at the switching instants ¢; and t; + At;. Since the Fourier series 
tends to the mean value of the function at a discontinuity at time, ¢, then a term 0.54 yi(t*) — Vi (¢} lim;—.0sG(s) with 


appropriate sign must be included in equation (6.4) to obtain the correct value of c;(t) at these two time instants. It is 


convenient to express c(f) and c(f) in terms of the A locus [6.6, 6.7], A(@,w), where 0 = wt, and is defined for a transfer 


function G(s) by 

Ac (0,@) = Re Ac(0,) + jlmAc(0,0) (6.5) 
with 

Re Ac(0,) = EVO, nw) sinn@ + Uc (0,nw) cos nO (6.6) 
and 

ImAc(0,@) = = (1/n)| Ve (0, nw) cos nO + Uc(0,nw) sin nO | (6.7) 


Ac(6,) isa generalised summed frequency locus with its real and imaginary values at a particular frequency, w, depending 
on weighted values, according to the choice of, 0, of the real and imaginary values of G(jw) at the frequencies nw for 
MND since co, In particular for @ = 0 the real (imaginary) part of Ac(@,@) depends only upon the real (imaginary) 
part of G(jw) at frequencies nw. It is possible to evaluate Re Ac(9,@) and Im Ac(@,@) in closed form for specific transfer 
functions from known infinite series, as is discussed in the next section, or approximately by computationally summing 


to a finite number, n, of terms. 
Using the above relationships in equations (6.2) and (6.4) gives 

ci(t) = (hAAtG(0)/T) + (h/z)[ImAc(— wt + wt,w) — ImAc(— wt + wt: + wAt,o) | (6.8) 
and 

c(t) = (wh/z)[ReAc(— wt + wt,@) — ReAc(- wt + wt; + wAt,o)| (6.9) 


Equation (6.9) may also be seen to be obtainable directly from equation (6.8) since from the definitions of equations (6.6) 
and (6.7) it can be seen that 


dim[Ac(,@)] _ 


10 — ReAc(6,a) (6.10) 
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Any periodic pulse waveform can be constructed from a summation of waveforms y;(t) with different origins and pulse 


widths. For example, the symmetrical square wave, y(t), of period T is given by 

y(t) = Dara) (6.11) 
with y,(¢) switching positive at time “1, y(t), switching negative at time ; and with # = f, + 0.5T and Af, = Ah = 0.5T. 
Choosing the time origin 4 = 0, then the output c(t) becomes 

c(t) = (2h/z)[ImAc(— wt,w) — ImAc(— wt + 0.57,0) | (6.12) 
Since 0.5w7T = 7, it is easily shown that this gives 

c(t) = (4h/x)[Im Ae (— at, a) | (6.13) 
and, similarly, 

¢(t) = (4@h/z)[Re Ab (— at, a) | (6.14) 
Here the A loci Aé are again defined by equations (6.5) to (6.7) but with the summations including odd terms only, that 


ism=1,3,5....... co, Before proceeding to show how limit cycles can be found using these expressions it is appropriate to 


summarise some properties of the A loci and show how they can be obtained for specific transfer functions. 
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6.3 Properties and Evaluation of A loci 


Limit cycles in relay systems 


Some further properties of the A loci, in addition to equation (6.10), which are easily found from the definitions are:- 


(1) The A’ locus for @ = 0 is identical with the Tsypkin locus [6.3], A(w), if lims..sG(s) = 0, apart from a 


constant factor. The relationship is 


A(@) = (4h/7)A°(0,@) 


The A® locus can thus be regarded as a generalised Tsypkin locus. 


(2) The A locus satisfies the superposition property, that is if the linear plant 
G(s) = G(s) + Go(s), then 


Ac(@,@) = Ac,(6,@) + Ac.(6,@) 


(3) If Gi(s) = G(s)e “, then 


Ac (@,@) = Ac(@ + wT,) 


(4) The loci are periodic in 6, with period 27, that is 


Ac(0,@) = Ac(@ + 27z,w) 


(5) For A® the periodicity is odd, that is 


AG(0,@) =— AG(— 0,0) 


(6.15) 


(6.16) 


(6.17) 


(6.18) 


(6.19) 


Since any plant transfer function can be written in terms of a summation of transfer functions having a real or complex 


pair of poles, use of equation (6.16) allows A loci for most transfer functions to be obtained in terms of the A loci of the 


few basic transfer functions given in Table 6.1 in the Appendix 6.10. In addition A loci for transfer functions with time 


delays can be obtained using equation (6.17). Expressions for the A loci of the transfer functions of Table 6.1 are given 


in Tables 6.2 and 6.3 respectively, for the A° and A loci. For the real pole situations the loci are expressible in terms of 


well known sine and cosine series. 
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6.4 Solving for Limit Cycles 


The results of the previous section have shown how a closed form expression can be obtained for the output of a linear 
transfer function with a periodic pulse waveform input. Thus the signal fed back to the relay input can be found. If this 


pulse type limit cycle is to exist then the following three conditions need to be satisfied:- 


(1) The input must pass through the relay switching values in the correct direction at the appropriate times to 
provide the assumed periodic output. The number of these ‘switching conditions’ to be satisfied will depend 
on the form of periodic mode being investigated. The simplest case is discussed in the next section for a 
relay with no dead zone with a symmetrical square wave output, which yields one switching condition and 
thus the requirement to solve one nonlinear algebraic equation. For complex waveforms there may be many 


switching conditions resulting in the requirement to solve many nonlinear algebraic equations 


(2) Once the switching conditions have been found the relay input waveform can be computed. This must satisfy 
the ‘continuity conditions, that is it must not pass through the relay switching levels at times which would 


cause the relay output to be different from that assumed. 
(3) The limit cycle must be stable. 


Although analytical solutions can be obtained for some simple cases of one or two switchings per limit cycle period, the 
best approach is to use software with graphics [6.8]. This then allows the above three conditions to be accomplished, 
namely solution of the nonlinear algebraic equations of (1), display of solution waveforms to check (2), and the stability 


computations, discussed in section 6.5, to be carried out for (3) 


6.4.1 Relay with no Dead Zone 


For a relay with no dead zone the simplest assumption for the periodic output is a square wave with 1:1 mark space ratio. 
Equations (6.13) and (6.14) give expressions for the waveform and its derivative produced at the output of the linear 
transfer function G, where the positive switching of the square wave was assumed to be at t = 0. Taking the input to 
the relay, x(t) =— c(t), then the required switching condition, assuming the relay to have a hysteresis of +A, that is it 


switches from -h to +h for an input of A and from +h to -h for an input of -A, is 
x(0°) = Aand x(07) > 0 (6.20) 


The ~ is used to show that the time should be taken before any possible discontinuity which may be created due to the 


switching. Thus using (6.13) and (6.14) one obtains 
ImAég(0,w) =—zA/4h and Re Aé(0,@) < 0 (6.21) 


provided lim,..sG(s) = 0. 
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For a given G(s) the solution for the unknown frequency of the limit cycle, w, can be found by solving the equality of 
equation (6.21) and checking the inequality, or done graphically from a plot of A¢(0,@) on a Nyquist diagram to obtain 
where it intersects the line — zA/4h. This approach allows easy comparison with the DF method where the solution is 


the intersection of G(jw) with the same line. Note from the definition A¢(0,@) has real and imaginary parts given by 


ReAi(0,@) = 7 Uc(0,nw) = S7ReG(jnw) (6.22) 
n=1(2) n=1(2) 
and 
ImAg(0,o) = >> (/n)V¥.(0,n@) = 9. (1/n)ImG(jno). (6.23) 
n=1(2) n=1(2) 


Thus the locus can be obtained approximately from these expressions by taking a small value for n, and exactly by summing 


to a large value of n computationally or using the analytical results of Table 6.2. 
For the specific example of G(s) = K/[s(1 + st)], one has using Table 6.2 that 

ImAé(0,@) = (Kzt/4)|(— 2/2A) + tanh(z/2A) | (6.24) 
where A = wrt. Thus using equation (6.21) the solution for the limit cycle is 

(x/2A) — tanh(z/2A) = AhK/t (6.25) 
and the inequality is satisfied as Re Ac(0,w) < 0 for all w. 
The DF method for this problem yields the equation 

A(1 +A’) = 4hKt/A (6.26) 
for the normalised frequency, A, of the limit cycle. 
For the particular case of K = 1,7 = 1, h/A = 3, the exact limit cycle frequency solution is 1.365 and that given by the 
DF method is 1.352. As expected if the solution is examined as the parameter, T, changes then, since as T reduces the 
transfer function is a worse low pass filter, the DF solution becomes less accurate but it is less than 3% in error for T > 0.25. 
When the square wave is asymmetrical, that is not a 1:1 mark to space ratio, the two pulse trains of equation (6.11) will 
be y.(t) switching positive at , = 0 with At, = Af, and y,(¢) switching negative at time f, = At and with At, = T — At. 


The required relay inputs of A and - A at the switching times of 0° (positive) and Ar (negative) then yield two nonlinear 


algebraic equations to solve for Ar and a. 
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6.4.2 Relay with Dead Zone 


To evaluate any odd symmetrical limit cycle in a loop having a relay with dead zone, the relay output can again be assumed 
to be the sum of two pulse trains as given in equation (6.11) but now with pulse widths Ar, = An = At ¢ 0.5T. The 
positive switching condition at time f,, which will be taken as 0, will be 

x(0°-)=0+Aand x(0-) > 0 (6.27) 
and the negative one at time Ar 


x(Atr) = 6+ Aand x(Ar) > 0 (6.28) 


The expressions for c(t) and ¢(f) can be shown to be 


c(t) = (2h/r)[ImAé (0, a) — ImA2(wAt, a) | (6.29) 
and 
C(t) = Qah/z)[Re AZ(0,w) — Re Aé(wAt,a)]. (6.30) 
3 
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Thus that the switching conditions, for x(t) =— c(t), are 


ImAé(0,w) — ImAé(wAt,w) = — (2/2h)(6 — A) 
ReAé(0,w) — ReAZ(wAt,w) < 0 


(6.31) 


and 


ImAé(0,w) — ImAé(wAt,w) = — (z/2h)(6 — A) 
ReAé(0,w) — ReAg(wAt,w) > 0 


(6.32) 


For a specific example one thus has two nonlinear algebraic equations (6.31) and (6.32) to solve for the two unknowns 
w and At. Once these are known the pulse waveform is known and the limit cycle waveform can be obtained from the 
A locus expression or by other means. Plotting this waveform then allows condition (3) in section (6.4) to be checked. 
To investigate the possibility of an asymmetrical limit cycle for the relay with dead zone requires the solution of four 


nonlinear algebraic equations. 


6.5 Limit Cycle Stability 


When discussing limit cycle stability using the DF method the perturbation equation (4.5), which has periodically time 


varying coefficients, was found and is repeated below. 
q(D) x(t) + p(D)n' (x * (t)) dx(t)) = 0 (6.33) 


In contrast to the DF case the exact limit cycle solution not a sinusoidal approximation is given by the above analysis for 


feedback loops with relays. Further if this equation is written in the state space form 

Ox * (t) = A(t) dx * (1) (6.34) 
then it can be shown that the periodically time varying matrix, A(t), is given by 

A(t) =A+ Ba'(x*(D)C (6.35) 
where (A, B, C) is a state space description of G(s). The matrix of equation (6.35) is piecewise constant and equal to 
A except over the infinitesimal time intervals when the relay switches. The value over the relay switching time may be 


calculated by taking the relay as a limiting case of saturation. 


For example if the relay switches from —h to +h this is the limit of an ideal saturation characteristic as € — 0 moving from 


a level -h at -e to +h at +e and the value of the second term in equation (6.35) is 


lim..oBC(h/e) = lim. .oBC[2h/x(t) A] (6.36) 
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where t, is the time at which switching takes place and A the duration of the switching time which tends to zero as ¢ tends 


to zero. Thus the product A(t)A over the infinitesimal switching time is given by 
A(t)A = BC[2hA/x(t)| (6.37) 


The term in the square brackets is the change in relay output divided by the derivative of the relay input at the switching 
instant. The value of this result is that a necessary and sufficient result, for the stability of a differential equation with a 
periodically time varying A matrix which is piecewise constant, is known[6.9, 6.10]. The result is that the system is stable if 
all the eigenvalues of the matrix Q have magnitude less than unity apart from one, corresponding to the periodic solution, 


which will have unit magnitude. The matrix, Q, is given by 
Q= [4.40 (6.38) 
int 
where the A; are the constant values of the A matrix over the time periods Az; and 
T=] [At 
i= 1 


is the limit cycle period with m values of the A matrix. For a symmetrical odd limit cycle Q may be computed over 0.5T. 
Thus for a limit cycle in a feedback system with an ideal relay with levels th switching at times 0 and 0.5T, the Q matrix 
is given by 

Q = exp[AT/2 ]exp[2hBC/x(0) | (6.39) 


6.6 Some Interesting Limit Cycle Problems 


For most feedback systems with relay elements the limit cycles obtained are of what might be called the ‘fundamental’ 
type discussed in section 6.4. That is (1) for a relay with no dead zone the relay will switch twice per period of the limit 
cycle, namely at times 0 and At in the period, T, where if the limit cycle has odd symmetry At = 0.5T and (2) for a 
relay with dead zone the relay will switch four times per period. In this section, rather than consider examples for these 
fundamental cases some more interesting examples are considered. The results for these problems were found using the 
theory described above and some extensions. The solutions were obtained using Fortran programs with interactive graphics 


to allow the display of the limit cycle waveforms. 


6.6.1 Example 1 - Invalid Continuity Condition 


Consider an ideal relay with outputs +1 in a negative feedback loop with the transfer function 
G(s) = (1 + 2s)/(s? + 1)(s + 1), [6.10-6.13]. The Nyquist locus of this transfer function starts in the first quadrant 
and moves out to infinity at unit frequency due to the undamped second order term in the denominator. After a phase 
change of -180° around an infinite semicircle it returns to the origin in the third quadrant. There is no intersection of the 
negative real axis so the describing function indicates no limit cycle. If the Tsypkin method is used then it can be shown 
that it gives limit cycle solution frequencies of m = 1/(2n + 0.5) for n = 1,2,.....cc. Checking the solution waveforms for 
these frequencies they are found to cross the zero switching level of the relay at times other than those assumed, due to 
the resonant nature of the plant transfer function. Thus they fail the continuity condition (2) given in section 6.4 and as 


predicted by the DF no limit cycle exists. 
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6.6.2 Example 2 - Multipulse Limit Cycles 


Here a very interesting practical example concerning the possibility of multipulse oscillations in the attitude control system 
of a communication satellite in a geostationary orbit is discussed [6.14-6.17]. The roll-yaw loop is controlled by a pseudo- 
rate controller which consists of jet thrusters with feedback through a time constant which has a different value when the 
thrusters are on to when they are off. Thus on its own the controller is a relay feedback loop containing a mode dependent 
transfer function. The satellite has solar panels for power collection which means that the plant transfer function consists 
of the rigid body dynamics of a double integrator plus lightly damped second order transfer functions, corresponding to 


the panel's resonant modes, in parallel. A block diagram for the system is shown in Figure 6.2. 
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Figure 6.2 Block diagram for system of example 2 (copied from reference 6.15) 
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Reference 6.15 discusses the determination of the multipulse limit cycles using the DF, but here the comments will refer 
to the use of the Tsypkin method discussed in papers [6.15, 6.16]. There are two original aspects in these papers the first 
being the extension of the Tsypkin method to evaluate quite complex limit cycles with multiple pulses per period and 
the second being an extension of the method to find limit cycles in mode dependent relay systems. Here the term mode 
dependent is used to refer to the fact that the loop transfer function changes with the relay output level. This topic is 


considered further in reference 6.17. 


To determine the possible existence of an odd symmetrical multipulse oscillation of n pulses per half period requires a 
solution to be found to 2n algebraic equations. Results were obtained by this approach in reference 6.16 for one assumed 
flexible mode and are given for the pseudo rate controller with equal and, then the usual case, unequal feedback time 


constants. Results with up to five pulses per half period were calculated and checked with analogue simulation results. 


6.6.3 Example 3 — Limit Cycle with a Sliding Mode 


The investigation of the step response of a relay control system in section 2.3.2 showed that a sliding motion, where the 
relay switches rapidly, was found. This therefore raises the interesting question of whether a limit cycle can exist with a 
sliding mode and if so how can it be determined theoretically. This problem is discussed in references 6.10, 6.18 and 6.19 
where a relay with a dead zone of +1, output levels +1 and linear transfer function G(s) = 8(s — 0.5)/(s° — 0.55 + 1) 
was studied. The relative degree of the transfer function was taken as one since it can be shown that sliding motion under 
relay control is not possible for a relative degree greater than unity. The Tsypkin method gave the solution for a limit 
cycle shown in Figure 6.3, which violates the continuity condition. The interesting point is that the waveform reverses 
immediately after the switching at A and then soon reverses again to pass through the switching level at C, which is argued 


is an indication of sliding with the first estimate for the time of sliding being from A to C. 
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Figure 6.3 Tsypkin solution for limit cycle (taken from reference 6.19) 
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To see if a sliding solution could be obtained the relay output was approximated by a straight line over this period and the 
Tsypkin method extended to obtain a solution for a relay output waveform having this shape. A new limit cycle solution 
was then found showing that after the switching the calculated relay input remained almost constant as required to verify 
sliding. By adding a further linear segment to approximate the average relay output waveform in the sliding mode by 
two linear segments the method gave the limit cycle shown in Figure 6.4 with a frequency of 0.386 rads/s, which was 
in agreement with the simulation result. It can be seen that the relay input is now almost constant at the switching level 


during the sliding period. 


Sliding region approximated by two Linear segments. 


Relay .nput 
xie! 


34 


8 
THETAC radians > 


-.34 


Figure 6.4 Limit cycle solution using Tsypkin extension (taken from reference 6.19) 


The example has also been done in simulink using a fixed step integration algorithm since the default option with a 
variable step length stops when sliding starts. The simulink diagram is shown in Figure 6.5, where the blocks in the upper 


left hand corner generate a time vector, z, of equal length to the output vectors x and y. 


Constant Integrator To Workspace1 


To Workspace 


Integrator] Integrator2} Gain 


Relay1 


Figure 6.5 Simulink block diagram for system having limit cycle with sliding 


The limit cycle at the input to the relay, x, is shown in Figure 6.6 and in the x-y plane in Figure 6.7. The initial condition 


was set slightly off the limit cycle. 
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Figure 6.6. Waveform of limit cycle at the input to the relay. 


Figure 6.7 Limit cycle plot in the phase plane 
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6.6.4 Example 4 - Multiple Limit Cycle Solutions 


For a feedback loop containing a single valued nonlinearity, which has a C(a) locus on the negative real axes and a linear 
transfer function which has a frequency response with more than one intersection of this axis, various possible limit cycles 
may result. A study of this situation using the approximate DF analysis, where both cubic and ideal relay nonlinearities are 
taken, can be found in reference 6.20. The Tsypkin method enables the existence of ‘fundamental’ limit cycles and their 
stability to be found but not limit cycles at combinations of the intersection frequencies. To illustrate the situation the 
case of an ideal relay in a negative feedback loop with the linear transfer function G(s) = 12(s + 1)7/s*(s* + 0.3As + A’) 
as considered in reference 6.20 is discussed. The Nyquist locus for the transfer function with \ = 4.0 is shown in Figure 


decreases as 


6.8. The negative real axis is cut at the two frequencies w, and w, and the ratio a =|G(jw) |/|G (jan) 
decreases. Using the IDF approach of the DF method the limit cycle at frequency w, is always unstable and that at a, 
is stable if a > 2 which corresponds to \ > 4.05. This compares with the Tsypkin solution which gave an unstable limit 
cycle at frequency w, for \ > 3.73, but no solution for smaller values of \. The limit cycle at w, was only stable for \ > 
3.73. Simulations gave good agreement with the Tsypkin theory for a single mode limit cycle but it cannot predict the 
combined mode. Using the SSDF and the ISSDF, the incremental describing function in the presence of two sinusoids, a 
stable combined frequency limit cycle was found to exist for \ < 4.05 in reference 6.20. Good agreement was also found 


with the amplitude and frequency components in the combined mode. 


Nyquist Diagram, 


Imaginary Axis 
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Figure 6.8 Nyquist plot for example 4 


6.6.5 Example 5 - Chaotic Motion 


Although no theoretical solution can be found for a chaotic motion there has been significant research done in trying to 
find conditions for its existence in autonomous feedback loops. In these studies both the approximate describing function 
approach [6.21] and Tsypkin method [6.22] have been used. Useful starting points are provided by the conjectures of 


Ogorzalek [6.23] and from the results of reference 6.21. Ogorzalek conjectured that for a negative feedback system with 


G(s) =— BI (s" + ni15 HF asssseese + Ms + Q) (6.40) 
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where all the coefficients are real and positive, and with nonlinearity n(x), chaotic motion was likely to exist if: 


(i) 


(ii) 


(iii) 


The Hurwitz sector of the linear system with n(x) replaced by K satisfies K, < K < K,, where K, = a/ — B 
and K; is finite. 


The nonlinearity n(x) intersects the line through the origin of slope K; at least once to ensure the 


existence of an equilibrium point. 


The slope of n(x) at all the intersection points with the lines of slope K, and K, should be either greater 


than the upper bound K, or less than the lower bound K>2, thus ensuring the instability of the equilibria. 


From their studies using the describing function Genesio and Tesi [6.21] concluded that the negative feedback system 


with nonlinearity n(x) and linear dynamics G(s) might be expected to exhibit chaos if 


(iv) 


(v) 


(vi) 


(vii) 


The application of the DF method to the system gives a stable predicted limit cycle x*(t). 


There exists an equilibrium point E, different from the one that generates the predicted limit cycle, 


which has one eigenvalue in the open right half plane and the others in the left half plane. 


The abscissa x, belongs to the range of the predicted limit cycle, that is, x*(t) = x, for some time f. 


As the filtering properties of G(s) deteriorate the motion may be expected to change from a limit cycle 


to chaotic motion. 


The investigations in [6.22] considered, n(x), to be a relay and used the transfer function form of equation (6.40) for G(s). 


It was conjectured that for n > 2 the system was likely to possess a region of chaotic motion if 


(i) 


(ii) 


Equation (6.40) has at least one pair of complex conjugate roots with positive real part 


The Hurwitz sector lies within the nonlinearity sector 


Based on these assumptions 3" and 4" order relay systems with chaos were synthesized. For the third order system with an 
p 'y SY: y’ yy’ 


ideal relay with unit levels and 6 = 1, then the upper and lower bounds of the Hurwitz sector are K; = a) and Ky = Q — QQ). 


Satisfaction of condition (iii) above requires 0 < K, < Ki, which gives a > @:Q@2. Thus if the transfer function is chosen 


with a@ = 4and a = 1.25 chaotic motion may exist for 0 < @ < 3.2. 
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The equilibria are stable for @ > 3.2, become neutrally stable with roots at -3.2 and +1.118j for a, = 3.2, and as q 
decreases the complex roots have positive real parts. Using the program based on the Tsypkin method no limit cycle was 
found for @ > 3.2 but several unstable limit cycle solutions were found as @, was decreased below 3.2. In particular an 
almost sinusoidal unstable limit cycle was found together with several ‘spiral’ type unstable limit cycles. When @ was 
reduced to 2.8 one of the ‘spiral’ limit cycles became stable and is shown in Figure 6.9. This may be referred to as a 3 spiral 
because of the number of oscillations. For 2.8 < @ < 3.2 the various unstable limit cycles had different numbers of spirals 
apart from the almost sinusoidal one which bounded all the others. Simulations starting with initial conditions inside 
this bounding limit cycle resulted in chaotic motion with the motion ‘jumping’ between different ‘spiral’ type unstable 


limit cycles as shown in Figure 6.10. 


Ase oie 
$9+2.852+) 25544. 
Nek TYPES 4 GAIN=-1 .0000 
N-L_ PARAMETERS TIME DELAT#0.0000 
1.0000 BIAS TO N-L 17P#0.0000 


RELAY INPUT WAVEFORM L-C Freqe0.19957 


Thete (rediens) 


Figure 6.9 Three spiral limit cycle (taken from reference 6.19) 
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Figure 6.10 Chaotic spiral type motion (taken from reference 19) 
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6.6.6 Example 6 - Unstable Limit Cycles and Simulation 


An interesting aspect of the above theory for finding limit cycles in relay systems is that both their waveform and stability 
can be accurately calculated. Stable limit cycles are normally not difficult to find in simulation studies as they may be 
reached from various initial conditions, unstable limit cycles are, however, much more difficult to locate. Since knowledge 
of the exact limit cycle solution allows initially conditions to be used in a simulation which lie on an unstable limit cycle, 
reference 6.23 studied the ability of a simulation to track the limit cycle. It was found that longer tracking was achieved, 
as expected, by using a smaller step size and also if the magnitude of the largest eigenvalue of the stability matrix Q was 


not significantly greater than unity. 


6.7. Forced oscillations 


The only modification required to be made to the method in order to investigate forced oscillations is that the relay input 
becomes a combination of the feedback signal and the forcing signal. Thus, if the relay input, x(f), is the error signal of 
the feedback loop then x(t) = r(f) — c(t). If a relay output waveform is assumed which switches when the input has a 
value of A at the switching time Af then r(At) — c(At) = A and x(¢) must have the appropriate sign. An example is given 
in reference 6.12 for the theoretical solutions to applying a sinusoidal input to a feedback system having a relay with dead 
zone. The system gain is such that it has a limit cycle with no input. A sinusoidal input with a frequency of roughly three 
times that of the limit cycle is applied to the system and its amplitude slowly increased. Assuming initially a relay output 
with switchings appropriate to the limit cycle the plot of the relay input as the input amplitude increases contains more 
of the input frequency which eventually leads to a solution with false switchings. It is then assumed that the relay output 
will have two rather than one pulse per half cycle and a solution confirmed. Thus, embedding the theory in software with 
graphical output enables some results to be obtained for forced responses both for systems which may be stable or possess 


a limit cycle in the autonomous state. 


6.8 Conclusions 


This chapter has shown how an analysis based on a method originally formulated by Tsypkin can be used to find the exact 
form of limit cycles in relay systems. In addition to evaluating the limit cycles the method also allows their stability to be 
found. The evaluation of the limit cycles, however, requires investigation of the solution of nonlinear algebraic equations 
whose number increases with the number of pulses per period at the output of the relay. This is not a problem if the 
approach is imbedded in computer programs with interactive graphics as described above for some interesting examples 
studied using software written in FORTRAN. It is also possible to extend the method to cover systems with more than 


one relay as well as relay feedback systems with harmonic inputs. 
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6.10 Appendix 


The appendix presents formula from which the A loci of various transfer functions can be evaluated. This is done by 
putting any transfer function into partial fractions in terms of the simple transfer functions given in Table 6.1 and then 


using the results for their A loci given in Tables 6.2 and 6.3. 


Type No. Transfer Function Type No. Transfer Function 
2 Wd+sN 


3 1/sT(. + sT) 8 1/11 +sT)' 
9 Is? + 2kseo, + @) 
5 1/(1 + sT)’ 10 wo! (s° + 2€sw, + >) 


Table 6.1 Basic Transfer Function Types. 


The results for transfer functions with real poles are defined in terms of sine and cosine series and for the complex poles 
in terms of some parameters. The results are given for the A® loci in Table 6.2. The definitions of the sine and cosine series 
and the results for the transfer functions are given in Tables 6.2(a) and 6.2(b). The parameters for evaluating the complex 


transfer functions are given in Table 6.2(c). 


Similar results for the A loci are given in Tables 6.3, 6.3(a), 6.3(b) and 6.3(c), respectively. 
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Table 6.2 Results for A° loci 


Download free ebooks at bookboon.com 


117 


An Introduction to Nonlinearity in Control Systems Limit cycles in relay systems 


4 
5° a g° (6, A) = n~ sin n8 
j,k j,k 22102) 4,2 2) 


for 0 <6<T 
SS ree -6,A) = “Sy 68d) 
$71,069) = w/4 


S®° (6,4) = mncosh[n-26)/2A]/4A* cosh(1/2A) 


° ae os re, 
B71 1899). = Sly 9 SA “y x 
268s A) = {16 sinh[(1-6)/2A] 
+ = sinh(6/i) /2 cosh(n/2A) W8A> 
x cosh(1/2i) 


- r75° 


° i o 
PA ee Le 


4 Fe 


S730) = {w"n sinh(6/A) + 278A cosh(1/2A) 
x sinh[(7-26)/2A] - on-6 cosh(8/i) 
+ 2107 cosh(n/2d) cosh{ (1-29) /2A] 
- n°sinh(6/A) tanh(s/2d) 1/64\°cosh=(4/2A) 


53,0 699%) = in“6 - 167) /8 


Table 6.2(a) Results for Sine Series 
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Table 6.2(b) Results for Cosine Series 
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Table 6.2(c) Parameters for Complex Pole Transfer Functions 
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Limit cycles in relay systems 
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Table 6.3 Results for A loci 
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Table 6.3(a) Results for Sine Series 
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Table 6.3(c) Parameters for Complex Pole Transfer Functions 


Download free ebooks at bookboon.com 


124 


Please click the advert 


An Introduction to Nonlinearity in Control Systems Controller Tuning from Relay Produced Limit Cycles 


7 Controller Tuning from Relay 
Produced Limit Cycles 


7.1. Introduction 


Control theory is typically taught from the starting point of having a mathematical model of a plant, most often linear, 
and a controller is then designed to provide a specified performance. Techniques usually regarded as classical control 
consider simple forms of controller, such as phase lead or PID and the problem then becomes, once the specific form 
of controller has been chosen, to try and find parameters for the controller so that the system performance meets the 
specifications. When working with a mathematical model of the plant many methods of classical control can be used to 
find the controller parameters. Dependent on the specifications there may be no solution, a unique solution, or many 


solutions to the problem. 


The latter, for example, could be the case if the only specification was on the response to a step input and required an 
overshoot less than 10% and a relatively slow settling time. On the other hand if it was required to minimise the integral 
squared error for the step input then there would bea unique or, possibly no feasible, solution. The most common adjustable 
parameter controller is the PID, with the three terms P- proportional, I - integral and D - derivative. It has been used in 
the process industry since the 1940’s first in pneumatic, then successively in vacuum tube, transistor, integrated circuit, 


and now usually in microprocessor form; an implementation that allows significantly more functionality to be included. 
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When commissioning or examining control loops engineers or plant operators often found that with suggested or existing 
PID parameter settings the loop performance was not satisfactory. Thus the obvious question arose, what could one do 
‘on the job’ to get better parameter settings. Enter Ziegler and Nichols [7.1] who in their 1948 paper suggested two, in 


principle, simple approaches:- 


1. Perform an open loop step response on the plant. From this response estimate the parameters of a first order 
plus dead time (FOPDT) model. Then set the controller parameters using given formulae involving the 


estimated plant parameters 


2. With the plant under closed loop control put the controller into the P mode. Adjust the gain until an 
oscillation results, measure its frequency and record the P gain. Then set the controller parameters using 


given formulae involving the oscillation frequency and the P gain. 


They referred to their approaches as parameter tuning, which seems logical as what was done, possibly by a technician, 
was perform a system test and from the readings calculate the controller parameters according to the given formulae, quite 
a standard engineering approach. Today many papers refer to parameter tuning when the values are calculated starting 
from a mathematical model of the plant, which does not seem the correct terminology, as one never sees the term used 


for obtaining the parameters of other fixed controllers, such as lead or lag controllers, by this approach. 


The Z-N methods consist of two parts, which is not often clearly stated, namely a testing ‘identification procedure and 


then the use of their tabulated formula to set the parameters. The ‘identificatiom procedures in the tests are:- 


1. An estimated FOPDT model from the step response. This was chosen because it was judged to be a 
reasonable model for many simple processes, although other models could be fitted to the response if 
felt more appropriate. It should be noted that in the early days of control the standard experimental 


identification tests were step and/or frequency responses. 


2. An estimate of the ‘so-called’ critical point (or gain margin point) of the plant, namely the gain and 
frequency value, w., when the plant gives a phase change of -180". This, is of course, a ‘small bit of 
information’ about the plant, and the major contribution of Z-N was to recognise the value of this 


information. 


Their suggested controller parameters based on the measurements from both tests were based on designs to give a relatively 


high overshoot, around 25%, in the closed loop step response. 


Given the identification result of 1, namely an FOPDT plant model, then there are basically an infinite number of analytical 
approaches which can be used to calculate PID controller parameters. The point is that the approach used should meet 
specified performance criteria which are application dependent and may, of course involve considerations related to input 
step and disturbance responses, noise, nonlinearity and other factors not usually mentioned in simple approaches. If the 
specifications address a desired closed loop step response behaviour, as many do, then an open loop frequency domain 


approach is indirect and in general will require iteration to achieve the desired closed loop performance. 
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Both the Z-N identification methods have problems in practice. In the first test noise in the measurement affects the accuracy 
of getting the FOPDT model and it cannot, of course, be used for an open loop unstable plant. If used as intended, however, 
for plants in which one has reasonable confidence have essentially the assumed dynamic behaviour, and not any ‘black 
box’ model, then the identified parameters normal provide a satisfactory model for analytical design methods [7.2]. The 
second test is extremely difficult to implement in practice as many process plants have slow dynamics, so the reaction to 
a change in loop gain is slow. Further the principle is based on linear dynamics and actuators invariably have dead zones 
and indeed without saturation of some form dangerous amplitude oscillations can result from too large an increase in the 
P gain. It is, of course, a test which can be performed ideally in simulation using a modern digital simulation language, 


although it was not easy on an analogue computer simulation and can be, as mentioned, a ‘nightmare’ in practice. 


Recognising these practical difficulties with the test, Astrom and Hagglund [7.3] suggested replacing the P gain control by 
an ideal relay, which became simple to implement with the move to microprocessor controllers in the 80’s. The question 
then arises, having replaced the P control by a relay, as to how one uses the resulting limit cycle to obtain information 


about the plant. The earlier chapters provide the information for doing this. 


7.2 Knowledge from the Limit Cycle 


This can be obtained by basing the analysis either on use of the DF method or the exact method for limit cycles covered 
in the previous chapter. The original proposal [7.3], usually known as relay autotuning was based on the DF method and 
like the original Z-N work assumed the plant was unknown but typically having a transfer function form common for 
the process industries. Recently other approaches have been suggested where by assuming a specific model form for the 


plant its parameters are estimated from measurements on the limit cycle. 


7.2.1 Relay Autotuning 


For a negative feedback loop containing a plant with transfer function, G(s), and relay, use of DF analysis gives 

N(a)G(Gja) =- 1 (7.1) 
which can be written 

| N(a)||GG@)| = 1 and ArgN(a) + ArgG (jw) =— 180° (7.2) 
For an on-off relay of height +h and hysteresis +A then 

N(a) = ae where 6 = sin '(A/a) (7.3) 


Thus when the nonlinearity is an ideal relay and denoting the limit cycle frequency by w,, one has 


ArgG(j@.) =— 180° and |G(Ga.) 


= am/4h (7.4) 
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This shows that the DF analysis yields 


= am/l4h = 1/K. (7.5) 


®. = w, and |G(ja.) 


for the frequency and magnitude of the critical point, where K, is the critical gain. Due to the approximation of the DF 


method this result is approximate and is normally more accurate the nearer the limit cycle is to being sinusoidal. 


It should also be noted that a is not the amplitude of the limit cycle at the relay input but the amplitude of its fundamental. 
Most practical implementations measure the peak to peak limit cycle amplitude and take a to be one half of this. It is 
possible to include additional computation software in the controller to obtain a better estimate of the fundamental a but 
in view of the inaccuracies in measurement it is probably not justifiable in many instances.. On the other hand sometimes 
easy corrections are possible [7.4] for example if the limit cycle is almost triangular, as it is for an FOPDT plant with 
small time delay, then one can calculate the fundamental from the peak value. Also the theory gives a better percentage 


accuracy for the frequency than the amplitude [7.5]. 
If time is no object in performing limit cycle experiments then one can:- 


(1) normally get more accurate measurements by averaging over several cycles. 

(2) vary the hysteresis in the relay to get additional results 

(3) include known extra linear elements in the loop to get additional results 

(4) get improved results from the DF method by making the limit cycle nearer to a sinusoid with the inclusion 


of a known tuned filter. 
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Typically (2) and (3) can be used to get additional points [7.6] to that of the critical point on the plant frequency response 
and iterative tuning of a resonant filter [7.6] can achieve (4) by producing a pure sinusoidal limit cycle at the plant critical 
frequency, which then yields the exact value of the critical point apart from, of course, measurement errors. The tuned 
filter was introduced in [7.6] in a paper concerned with the design of nonlinear PID controllers for nonlinear plants, where 
to get better information on the plant frequency response a sinusoidal input was required. A big advantage of the relay 
autotuning method is that the amplitude of the limit cycle can be controlled by the setting of the relay output levels, and 
this avoids one of the practical problems of the second method of Z-N. Further one can use this feature to estimate whether 


the plant behaves reasonably linearly and if not use the results to design nonlinear PID terms as discussed in section 9.5. 


7.2.2 Plant Parameter Identification 


It has been shown how the limit cycle in a loop with a known plant transfer function under relay control can be calculated 
exactly using the Tsypkin method or estimated using the DF method. In particular for the case of an odd symmetrical limit 
cycle the first method allows exact determination of the limit cycle amplitude and frequency at the relay input in terms of 
the plant parameters. Thus, the method can also be used for plant transfer function identification since measurement of 
the amplitude and frequency of the limit cycle will allow two unknown plant parameters to be found from the equations, 
for example the time delay and time constant in an FOPDT model when the plant d.c. gain is known. Analysis for an 
asymmetrical limit cycle and measurement of more parameters in the resulting limit cycle enables additional plant 


parameters to be calculated, for example four parameters in an SOPDT model. 


Several papers [7.7, 7.8] describe the identification of plant parameters based on an assumed plant transfer function 
form, using exact limit cycle analysis and limit cycle measurements. From a practical viewpoint they have the following 


possible disadvantages. 


(1) The microprocessor controller software, as well as being able to ‘measure’ various limit cycle parameters, 
must be able to solve the appropriate nonlinear algebraic equations to get the plant parameters. Further it 


must contain a design procedure for calculating the controller parameters for the known plant. 


(2) There will be errors in the measured limit cycle parameters, due to noise etc, and these may cause problems 
in solving the nonlinear algebraic equations, too large errors for example can result in no or physically 


unallowable solutions. 


It is also possible to assume a plant transfer function and determine expressions for some limit cycle parameters using DF 
analysis rather than the exact approach. Again one can then use the equations ‘in reverse’ and calculate plant parameters 
from limit cycle measurements. If this is done then due to the approximation of the DF method, particularly when dealing 


with asymmetrical limit cycles, the problem cited in (2) is normally worse. 
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7.3 Tuning the Controller 


Many design methods may be used to obtain controller parameters to meet given system specifications when the plant 
transfer function has been estimated by one of the methods in section 7.2.2. Both ‘classical and ‘ modern methods’ will 
often turn out to require an iterative approach as they will invariably not address the required specifications directly 
but by inference. For example, if the design specification is with respect to the input step response, then if one uses a 
frequency domain design one will be relying on an imprecise link, say between phase margin and overshoot in the step 
response, or if one uses pole placement the effect of any closed loop zeros or neglected poles on the expected response 
will have to be checked. 


The actual system behaviour will differ from the theoretical calculations due to the modelling errors of the identification 
and this, assuming the plant behaves reasonably linearly, can be taken into account by considering parameter uncertainty. 
Much is made about the concepts of robust control for such situations but no results address the real design problem, 
namely if the plant parameters vary over a certain range will the specifications, which normally relate to responses, not 
stability, still be met. Fortunately today MATLAB is so fast that the robustness of a controller design can easily be checked 
simply by doing a set of Monte Carlo simulations over the expected range of variations in a few parameters. In fact a 
good designer may be able to identify the ‘worst case; so that the design just meets the specification for this case and is 


better for other parameter values. 


On the other hand if one uses the limit cycle test to estimate the critical point of the plant transfer function, two issues 
arise namely how accurate is the critical point estimate and secondly how, knowing just this bit of information, can 
the controller parameters be found so that the system meets the desired specifications. The first point has largely been 
addressed in the previous sections and to some extent how errors in the critical point value may affect the subsequent 
design can be examined by considering in simulation how the performance differs when selecting a few slightly different 


values for the critical point. 


The second point is much more difficult, since even if a simple step response specification is given, one cannot assert for 
given PID parameters that if a plant with a given critical point has a closed loop step response with a given percentage 
overshoot then all plants with the same critical point will have the same overshoot. Certainly the closeness of the open 
loop frequency response of the compensated system to the (-1, j0) point has a significant effect on the form of the 
closed loop step response but the shape of the response in this region can differ significantly for different plant model 
transfer functions with the same critical point and the same compensator. It is therefore desirable in using critical point 
information to have some idea of the form of the plant transfer function. Having obtained the critical point all that can 
be done in a controller design algorithm is to select the controller parameters so that, subject to realisability constraints, 


the compensated frequency response has the critical frequency at a desired point in the Nyquist plane. 


This point can be moved around, within a finite region, by adjusting the controller parameters. For the ideal PID controller 
K[1 + s7+ (1/sTi], Z-N chose to take 7; = 47; and it is easily shown that their choice for T, and K places the critical 
frequency on the compensated locus at a gain of 0.66 and phase angle of -155°. The corresponding point for their suggested 
PI parameters is 0.46 at a phase angle of -191°. These choices, as for their first method, tend to produce closed loop step 


responses with a relatively high overshoot. 
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Thus the success of PID controller design based on critical point information is one of knowing where the critical frequency 
should be on the compensated locus to, hopefully, meet the design specifications. Given some idea of the plant transfer 
function it is relatively straightforward to come up with good strategies. For example, in reference 7.9 optimum parameters 
are found for PID controllers to minimise integral of error performance indices for a step input to a closed loop with an 
FOPDT plant. These are given for minimization of the ISE, ISTE and IST’E. Typically as the time weighting increases the 
overshoot decreases, the rise time is slower and the settling time increases slightly. The location of the critical frequency 
on the compensated Nyquist locus has also been calculated for these designs so that critical frequency shifting designs 
can be used for other similar plant transfer functions. Minimization of the ITSE typically leads to a design with around 
10% overshoot. Many different design approaches are possible for deciding where to place the critical frequency. Also, 


depending on the application different requirement specifications may have to be met. 


One might be tempted to say that knowing a bit of information about a plant transfer function, namely its critical point, 
cannot be anything like as good as having an experimentally measured open loop plant step response, in principle a 
multiplicity of points. The problem is that measured step response results are not exact with large amounts of noise 
often present in industrial measured data. To illustrate a point consider the two step responses in Figure 7.1, curve (a) 
is for a chosen transfer function and curve (b) for a transfer function derived from it by a model reduction technique. 
On this basis, particularly if the response of (a) had resulted from measurements contaminated by noise, one would 
have clearly concluded that the transfer function of curve (b) would be a good model to use for controller design for the 
original transfer function (a). Now go to Figure 7.2 and look at the corresponding frequency responses. It can be seen 
that our previous conclusion is clearly not valid as the transfer function with the model (a) has a finite and model (b) an 
infinite gain margin. Thus a poor design for say, a constant gain controller, would result from using model (b), one in 


fact which would probably result in an unstable system with the chosen transfer function. The two transfer functions are 
13s + 10 aiid 5.435 
(s + 1)(0.5s + 1)(s? + 35 + 1) s+ 1.7395 + 0.5435 


, respectively. 


Figure 7.1 Step responses of original and reduced models. 
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Figure 7.2 Frequency responses of original and reduced models 
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7.3.1 An Example 


K ee 
(s + 1)? 


Consider for example a plant with the transfer function G(s) = with K, = 1. 

The objective is to show that good tuning can be obtained for a PID controller based on an autotuning test. Using a relay 
of unit height the measured values of the amplitude of the limit cycle, A, and frequency, @, from a simulation with 
no noise introduced were A = 0.287 and @ = 1.90. Using equation (7.5) with a = A gives w. = 1.90 and K. = 4.44 
which compare with the actual values of 1.92 rads/s and 4.69 for the known plant transfer function. The normalized gain 
«x = K,K. and assuming K, = 1 has been found by other measurements « = 4.44. The tuning formula, suggested in 
reference 8, which gives optimum parameters for ISTE tuning (minimization of the integral squared of time multiplied by 
error) of a FOPDT plant based on x,K. and w., gives K = 2.55, 7’ = 2.61 and T = 0.40. The unit step set point response 
for these tuning parameters in the controller is shown in Figure 7.3 (response (1)) together with that for controller 
parameters of K = 2.42, 7 = 1.99 and 7; = 0.63 obtained by minimization of the ISTE criterion (response (2)) for this 
known plant transfer function. The results are very similar considering the approximations involved and the simplicity 


of the autotuning approach. 


Figure 7.3 Step response comparison for autotuning (1) and ISTE minimization (2) 


7.4 Autotuning using the Relay in Parallel with the Controller 


The standard autotuning approach is to replace the controller by a relay so that the forward path of the feedback loop 
contains a relay and the plant. Further if one wishes to check certain aspects of the compensated system, for example find 
its gain margin, then this can be done by placing the relay in series with the controller and plant [7.10]. These approaches 
involve what might be called separate tuning experiments but in some cases it may be appropriate to try and tune the 
system whilst it is still operating normally. This can be done by placing the relay in parallel with the controller [7.11]. 
This procedure was first suggested in reference 7.12, although in the paper the approach was used for plant parameter 


identification rather than critical point tuning. 
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When the relay is placed in parallel with the controller the feedback loop is as shown in Figure 7.4. To determine a limit 
cycle by DF analysis one simply opens the loop at the relay and calculates the transfer function, -G,, between e,, the relay 


output, and e, the relay input, in the figure. 


r(t) 


Figure 7.4 Relay in Parallel with the Controller 


The basic equations (7.1-7.2) still apply but with G replaced by G., It is easily shown that 


_ = G 
Gee 1+GG 70) 


which can be written 


fy CG. J 
G.14+GG. G 


G. = x Tr (7.7) 
where 7:, is the closed loop transfer function. Now consider the Nyquist diagram of GG,, and denote its magnitude at 
the frequency of the limit cycle, w,, by M and that of 1 + GG., by 0M, as shown in Figure 7.5. The magnitude equation 


in (7.2) can now be written 


1 M _ax 4h 


[GoM 7 4h which gives 0 = ae. (7.8) 
and the corresponding phase equation is 
b = 180° + Arg (1/G,) (7.9) 


where — ¢ is the closed loop phase angle shown in Figure 7.4. 
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Figure 7.5 Nyquist diagram and triangle for calculations. 


Since the transfer function in the controller is known |G| and Arg G_ are known at the frequency, w,, so that equations 
(7.8) and (7.9) can be solved, respectively, for p, as a is known from the limit cycle measurement, and ¢. Using the cosine 
rule with the angle ¢ for the triangle, since its sides are 1, M and pM, the value of M can be found. Arg(GG.) can also be 
found from the triangle so that both the magnitude and phase of the frequency response of G at the frequency , w,, can be 
calculated. This will only be the critical frequency for G if Arg G_ is zero, that is if the test can be done with the controller 
in the proportional mode. However, if w, is near to the critical frequency then settings based on the critical frequency 
should produce satisfactory results. Alternatively, other theories can be produced for finding controller parameters with 


various types of plant for other frequencies on the plant frequency response locus. 
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7.4.1 An Example 


— 0.58 
é 


As an example results from simulations of the FOPDT plant with transfer function G(s) = a4 


K(1 + 28h)’ 
AsT; 


and a relay with h = 1 
in parallel with the PID controller having 7; = 47), so that G.(s) = , are given. 

The critical frequency and gain of this plant are, respectively, 3.67 rads/s and 3.81 and for the simulation the parameters 
of the controller were set at the Z-N values of K = 2.284 and T, = 0.214. For the input step response this results in an 
overshoot of 52.5%. The amplitude and frequency of the limit cycle at the input to the controller, which is shown in Figure 
7.6, were measured as 0.680 and 4.30rads/s. 


SS 
O85 13.2 134 136 138 14 142 144 146 148 15 


Figure 7.6 Plot of limit cycle 


For the controller at w = 4.30, |G.(jw)| = 2.72 and ArgG.(jw) = 33.0° giving p = 0.688 and g = 147° from eqns. 
(7.8) and (7.9), respectively. Using the cosine formula for the triangle of Figure 7.5 gives M the magnitude of GG_ at 
w = 3.40 as 0.617, and applying the formula again gives a, the angle below the real axis of M equal to 13.4’, so that 
ArgGG_ = -166.6". Using the known values for the magnitude and argument of G,, then gives the magnitude and argument 
of G as 0.227 and -199.6°. Despite the fact that the waveform is not a good sinusoid these values are within 1% of the 
calculated values from the transfer function. The problem now becomes how should this information be used to set the 


controller parameters? 


Here use is made of the ISTE critical point design procedure of reference 7.8. For this plant the critical point needs to 
be moved to 0.58 at -153°, which means the PID gives a phase lead of 27°. Thus for a higher frequency it will give more 
phase lead, so that a reasonable suggestion for the frequency 3.40 is to increase the shift by around 25% to 35°. This gives 
a phase of approximately -165° (199.6° -35°) for this frequency on the compensated frequency response locus. Since this 
point is nearer (-1, 0) than the critical frequency point by 12° a slight reduction in the magnitude by around 10% to the 
value of 0.53 is taken. Using the same T/T, value as for the critical point shift design of reference 7.8, the required PID 
parameters are K = 1.94, T, = 0.205 and T, = 1.152, respectively. The input step response for this design is shown in Figure 
7.7 and as expected is compatible with an ISTE design. 
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Figure 7.7 Closed loop step response 


7.5 Conclusions 


This chapter has shown how information can be obtained about the linear dynamics of a feedback loop from measurements 
on limit cycles produced by adding a relay. In particular, the focus has been placed on using the results to set the parameters 
of a PID controller, however the concept can be used for deciding on appropriate parameters for other fixed structure 
controllers. With modern microprocessor controllers the approach is relatively easy and provides a good method for the 
automatic tuning of controller parameters without the need for mathematical modelling of the plant. Care, however, has 


to be taken to ensure that safe plant operation will be possible under relay control. 
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8 Absolute Stability Results 


8.1 Introduction 


It has been seen in previous chapters that the behaviour of the nonlinear feedback loop of Figure 3.1 and indeed, its stability, 
may be affected by the input. This has therefore led to two definitions of stability for a nonlinear system, the first stability 
of the autonomous system, that is with r(t) = 0 in Figure 3.1, which is often referred to as stability of equilibrium and to 
which most early work was directed. The second form of stability is known as bounded-input bounded-output (BIBO) 
stability and requires the output to be bounded for a bounded-input signal. These two forms of stability are, of course, 


the same for a linear system; for which necessary and sufficient conditions for stability have been known for many years. 


On the other hand the general nonlinear problem for the stability of equilibrium remains unsolved. Indeed until the 
development of frequency response based criteria in the 1960’s [8.1] the only available results were those from the work 
of Lyapunov which were not easily usable by an engineer. The frequency response based criteria, which will be presented, 
only give sufficient conditions for stability and thus often produce very conservative results for a specific nonlinear 
characteristic as they only use limited information on the nonlinearity. On the other hand a mathematician may argue 
that the results are very robust because they ensure stability for a large number of nonlinearities which satisfy certain 
conditions. This, of course, means that the system can be proven to remain stable if the conditions are not violated by 


any perturbations to the nonlinearity. 


The primary concern of this chapter is to cover frequency response absolute stability results as they allow comparisons 
to be made with both the DF method and limit cycle results for relay systems. In the next section Lyapunov’s method 
is briefly covered as some of the original frequency domain results were obtained using it whilst others were obtained 


using functional analysis. 


8.2 Lyapunov’s Method 


The direct or second method of Lyapunov for determining the stability of an autonomous system is a timedomain method 
using a state space description. The problem is to determine the stability of an equilibrium state (singular point), which 


without loss of generality may be assumed to be the origin of thestate space, for the system 
x = f(x) (8.1) 


Before considering Lyapunov’s method some preliminary topics must be considered. 
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8.2.1 Definitions of Stability 
Definition 1 


The equilibrium state 0, or the equilibrium solution x(t) = 0, is called stable if for any given positive e, there exists a 
positive 5, such that || x(t) || << 6 implies |]x(4) || < € for all t > t, where || || denotes the Euclidian norm. It should be 
noted that a system exhibiting a limit cycle or solutions for which|| x(¢) || increases temporarily will be considered stable 


by this definition ,when e and 6 are appropriately chosen, but it rules out solutions which grow without bound. For a 


nonlinear system stability may only exist in some domain || x || < R of the state space. 


Definition 2 


The equilibrium state 0 is called asymptotically stable if it is stable and x(t) — 0 as t +20 


Definition 3 


The equilibrium state 0 is called asymptotically stable in the large, or globally asymptotically stable, if it is asymptotically 


stable and the domain R covers the whole state space. 


What do you see? 


A skateboarder? 


We see an All Options Junior Trader leaving 
the office after a successful day. 
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8.2.2 Positive Definite Functions 
Let V(x) be a time invariant scalar function of the vector x and S be a closed bounded region in the x space containing 


the origin. 


Definition 4 


The function V(x) is positive definite (semi-definite) in S if, for all x in S 


(i) V(x) has continuous partial derivatives with respect to the components of the vector x 
(ii) V(0) = 0 
(iii) V(x) > 0 for x # 0 (or V(x) = 0). 


If the inequality signs in (iii) are reversed then the conditions for a negative definite (or semi-definite) function are obtained. 
Thus, for example, in two space if V(x) = x7 + x7 then it is positive definite as it is only zero at the origin, but the function 


V(x) = (%1 + x)’ is only positive semi-definite as it is zero everywhere along the line x + x = 0. 


A common quadratic form scalar function is 
Q(x) = x"Qx (8.2) 


where the matrix Q, which is nxn for x an n vector, is chosen to be symmetric. In this case it is known that a necessary 
and sufficient condition for the quadratic form of (8.2) to be positive definite is that the determinants of all the principle 
leading minors of Q be positive. For it to be positive semi-definite all the principle minors, not just the leading ones, 


must be non-negative. 


8.2.3 Lyapunov’s Theorems 
Theorem 1 


The null solution, or equilibrium state of the origin, of the system of equation (8.1) is (asymptotically) stable if in some 
neighbourhood of the origin there exists a positive definite function V(x) such that its derivative V(x) with respect to the 
solutions of equation (8.1) is negative (definite) semi-definite in that region. The conditions for asymptotic stability may 


be relaxed slightly by a theorem due to LaSalle which is given next. 


Theorem 2 


The null solution of the system of equation (8.1) is asymptotically stable if in some neighbourhood R ofthe origin there 
exists a positive definite function V(x) such that its derivative V(x) is negative semidefinitein R and is non-zero in R 
along any solution of equation (8.1) other than the null solution. If the conditions of Theorem 2 are satisfied for all x 
then the equilibrium state is asymptotically stable in the large. The function V(x), which satisfies the stable condition of 
Theorem 1, for the equation (8.1) is known as a Lyapunov function for the nonlinear equation. It should be noted that 
finding a Lyapunov function provides a sufficient condition for stability of equation (8.1). Thus the difficulty of applying 
Lyapunov’s method is one of finding a ‘good’ V(x). Failure to find one does not mean the system is unstable and finding 


one may result in a very conservative condition for stability. 
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8.3. Application of Lyapunov’s Method 


As might be expected it is possible to prove all the known results for the stability of a linear system using Lyapunov’s 


method but good results for nonlinear systems, despite intensive research by mathematicians, are not easy to obtain. 


8.3.1 Linear System 


For the linear system 

x = Ax (8.3) 
one can choose V(x) = x’ Qx. Differentiating and substituting for x and x” gives 

V(x) = x7(A7Q + QA)x = x" Px (8.4) 
where 

P=A'Q+QA (8.5) 
which is known as the Lyapunov equation, and which can be solved for Q given A and P. 


0 1 
A Lyapunov function for a linear system can be found by taking P = -I. For example if A = (! 9 3} then whit P = -I 


; : 5/4 1/4 
the solution for Qis Q = Ge ei, and 
V(x) = (1/4) (Sx? + 2x22 + x3) (8.6) 


which is positive definite. 


On the other hand if we had assumed any quadratic form for V(x) it may not have been possible to prove stability. For 


instance for 

v(x) = xp +.x3 (8.7) 
Then 
V(x) = Qin + Qxox = Qxix2 + 2x2(— 2x, — 3x2) = — (201% + 6x2) which is not negative definite and no conclusion can be 
drawn about stability. A physical interpretation can be found for this if the equations (8.6) and (8.7) are drawn on a phase 


plane for different constant values of V together with trajectories for the system. It will be seen that all the trajectories for 


the system cut the former constant V curves only once in the direction of decreasing V, which is not the case with the latter. 
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8.3.2 Nonlinear System 


There have been many efforts to try and find suitable V functions for nonlinear systems. One of the most successful was 


due to Lure who for a feedback system with a single nonlinearity in the feedback path given by 


x = Ax + bu 
u = n(z) (8.8) 
Z=cx 


used a V function of the form 
V(x,z) = x" Ox + B [ n@az (8.9) 


He was able to obtain necessary conditions for V(x, z) to be a Lyapunov function with (a) the nonlinearity confined to 
the sector [0, K] and all the eigenvalues of A in the left half plane and (b) the nonlinearity confined to the sector [¢, K] 
and all the eigenvalues of A confined to the left half plane or on the imaginary axis. Perhaps most importantly it led to 
results being found in the frequency domain which are more easily used by the engineer. 

8.4. Definitions and Loop Transformations 

Before presenting some frequency domain results some definitions need to be given andalso two useful loop transformations. 
The results have been derived from either or both Lyapunov’s method and functional analysis. 


Definition 5 


A system is bounded input - bounded output (BIBO) stable when a bounded input results in a bounded output. Typically 


a bounded signal is defined in terms of its L,-norm, that is 


lx = / fro dt , which is finite for a non-periodic BIBO signal. 
0 


Definition 6 


The A matrix of a transfer function G(s) is Hurwitz if all of its eigenvalues lie in the left hand side of the s-plane. This 
will be denoted A € {Ai} 


Definition 7 


If the A matrix is Hurwitz apart from one eigenvalue at s = 0 then A € {Ay} 


Definition 8 


If the closed loop system of Figure 3.1 has only eigenvalues in the left hand side s-plane , when the nonlinearity is replaced 
by a linear gain k, this will be denoted by Ac € {Ai} 
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Definition 9 - Sector bounded nonlinearity 

A nonlinearity n(x) which is a real, continuous, single valued, scalar function of x, with n(0) = 0 belongs to the class {N} 
if ki < n(x)/x Sk: for x £ 0 and is denoted by n(x) €[ki, kr]. 

Definition 10 - Slope bounded nonlinearity 

A nonlinearity n(x) which is a real, continuous, single valued, scalar function of x, with n(0) = 0 belongs to the class {M} 
if m < dn(x)/dx < m for x # 0 and is denoted by n’ (x) € [m,m | 

Definition 11 - Monotonic nonlinearity 

A nonlinearity belonging to class- {M} which also has m S[n(x1) — n(x2)]/[x1 — x2] S mp for all finite x,and x,#x, belongs 
the class {M, } 

8.4.2 Loop Transformations 


Two interesting transformations which are useful for obtaining stability results with the above criteria, and have also been 


used in the proofs of some of them, are the pole and zero transformations. 
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8.4.2.1 Pole transformation 


The closed loop properties for the autonomous system consisting of a nonlinearity n(x) and plant transfer function G(s) 
in the forward path are not changed if a negative feedback path with gain p is placed around the plant and a feedforward 


path of gain - p is placed around the nonlinearity as shown in Figure 8.1. 


Figure 8.1 Illustration of pole transformation for loop with no input. 


G(s) 
1+ 0G(s) 
Np(x) € [ki — 0,k2 — co]. The transformation as well as changing the poles of the linear dynamics changes the sector of 


The new autonomous loop has linear dynamics G,(s) = and if n(x)E[k,k| the new nonlinearity 


the nonlinearity but not the sector width defined by k,-k, for n(x) €[ki, kx] 


8.4.2.2 Zero transformation 


The added paths of gain o in the zero transformation are the reverse of those in the pole transformation, that is the 


feedback path is around the nonlinearity and the feedforward path in parallel with the linear dynamics. The new loop 
k, ky 
1+ 0k’ 1+ 0k 


elements are G,(s) and n,(x) where G,(s) = G(s) — o and n,(x) © 


8.5 Frequency Domain Criteria 


Several frequency domain criteria are given in the following sections. 


8.5.1 Popov Theorem 


A sufficient condition for the autonomous system of Figure 3.1 to be asymptotically stable with 


A. € {Ai} for k = 0 and n(x) €[0,k] is that there exists a non-negative q such that 
Re{(1 + jog)GGo)}+k'=d>0 (8.10) 


where 6 is an arbitrary small constant. When A; € {Ao}for k = 0 then the result is true for n(x) €[e,k] where ¢ is a 


small positive constant. Equation (8.10) can be written as 
Re(G(jw) > — k' + qaImG(jw) (8.11) 


If the real and imaginary axes of the Nyquist plot are denoted by x and y then this means that for stability the Nyquist 
plot of GGw) should lie to the right of the straight line 


x=—k'+quy 
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which is a function of the frequency w, so the line has to be changed for each frequency on G(jw) . The result is usually 


interpreted graphically in the frequency domain using the Popov locus 


G(jw) = ReG(ja) + jw) + jw ImG(jw) 


The theorem can then be stated as, ‘A sufficient condition for the system to be asymptotically stable is that the Popov locus 
G(jw) lies to the right of a straight line with finite slope passing through the point (-k',0), provided lim a) Pak 


The slope of the line is q' . 


The graphical interpretation of the Popov criterion is shown in Figure 8.2 for two different G(jw) loci labelled G1 and 
G2. Assuming the transfer functions are such that the linear system will be stable for gains, k, in the sector k €[0,k:] 
, which is known as the Hurwitz sector, then it is seen that the Hurwitz and Popov sectors are the same for G2 but the 


latter is smaller for G1. 


Circle for slope|sector 1-2 


Stable for both G1 and G2. 
Hurwitz and Popov sectors 
same for G2 but not G1 


Nyquist plot 


Figure 8.2 Illustration of the Popov criterion. 


8.5.2 The Circle Criterion 


There are two forms of a circle criterion which can be derived for a sector bounded nonlinearity from the Popov criterion. 
One is the case when q is taken equal to zero and the criterion is then usually referred to as ‘the circle criterion. For this 
case it has further been shown to be true for BIBO stability. It is shown in reference [8.1] that if a pole transformation 
is applied to the autonomous system, which results in a sector transformation for the nonlinearity and this is taken as 
n(x) €{N} then the Popov criterion for a stable G(s) is satisfied if the Nyquist plot of G(jw) does not enter the circle 
with centre at the point — l ( l + | 1 ge 1 _1 and which passes through the points (— k;',0) and (— kz', 0). The 


2\k  k,} 2 kk 
dependence of the circle on w is sometimes not a problem since q can be chosen large enough to take the circle centre 


well into the third quadrant for any relevant values of w. If q is taken equal to zero, however, which is a special and more 


restrictive case of the Popov criterion, then the result is the following theorem. 
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Theorem 1 - The Circle Criterion 


The system will be asymptotically stable for n(x) € {N} and a linear transfer function G(s), here the 


circle Dhas its centre on the negative real axis of the Nyquist diagram and passes through the points (— k;',0) and (— kz', O)if 
(i) for 0 < ki < k,the Nyquist plot of G(jw) encircles the circle D, P times counter-clockwise, where P is the number of 


poles of G(s) with positive real parts 


(ii) for 0 =< ki < k& the Nyquist plot of G(jw) satisfies ReG (jw) > — ky'and P = 0. 


(iii) for ki < 0 < k, the Nyquist plot of G(jw) lies in the interior of the circle, D. 


The criterion corresponds to the Nyquist criterion for linear systems for (i) with the Nyquist critical point (-1, 0) becoming 
the circle D. Figure 8.3 illustrates the situation for two stable transfer functions G1 and G2 whose Nyquist plots close 
with a semicircle in the right half plane and the circle is for a nonlinearity in the sector [1,3]. The closed loop nonlinear 
system with G1 is stable according to the criterion but no conclusion can be reached about G2. It should be noted that 
as k, tends to zero the circle tends to a vertical straight line through — k;'. Obviously, a disadvantage of the criterion for 


autonomous systems, since it is a special case of the Popov criterion, is that it produces more conservative results. 
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Circle, D, for sector 1-3 


A A Nyquist plot 


Figure 8.3 Illustration of the circle criterion 


It is interesting to compare this result with the DF method as for any nonlinearity n(x) € {NV} it was shown in section 
3.3 that V(a) € {N} so the DF for the class N of single valued nonlinearities is the diameter of the circle of Fig 8.3. Thus 
for ‘DF stability’ G(jw) should not cut the diameter of the circle. Aizermann conjectured that ‘DF stability’ was correct 
but the ‘so-called’ counter example of Fitts given in section 5. 5 shows this to be untrue. The feature of the Fitts example 
which breaks ‘DF stability’ down is that G(jw) exists in both the second and third quadrants (defined in a counter clockwise 
direction from the first quadrant within the positive axes). It is interesting to note that if G(jw) is confined only to the 
third and fourth quadrants then ‘DF stability’ can be proved by the more general circle criterion discussed at the start of 
this section. Since this is a typical region for many Nyquist loci, that is the plant transfer function has a phase lag between 


0 and 180°, this result may be regarded as the small phase shift stability theorem. 


8.5.3 The Off-axis Circle Criterion 

The off-axis circle criterion applies when the nonlinearity is monotonic. Its applicability is therefore restricted but it can 
give less conservative results than the other criteria. 

Theorem 2 


The autonomous system with G(s) such that A, and A,.2 € Aiand n(x) € {M,,}with n’ (x) € [1m 7] is asymptotically stable 
if the Nyquist locus G(jw) for @ €[0,co]does not enter the circle D> where D’ is a circle which passes through the points 


(— m;',0) and (— mz',0) with centre anywhere in the Nyquist plane. 


Remark 


For the case of m,=0 the circle becomes a straight line through (— m2',0), so that for stability the Nyquist locus must 
lie to the right of this line. Figure 8.4 illustrates the situation showing the loop stable for G1 but not for G2. However 


stability for G2 can be proved by choosing a larger radius for the circle and moving its centre farther away from the x axis. 
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Circle, D', for slope secto} 


Nyquist plot 


Figure 8.4 Illustration of the off-axis circle criterion. 


Remark 


Since the diameter of the circle can be made very large the limit of the chord is the line between the points (— m;',0) and 
(— mz',0) . Thus for a monotonic nonlinearity instability can only occur if the Nyquist locus exists in both the second 
and third quadrants. This again shows why a G(jw) like that in the Fitts counterexample breaks down a DF stability type 
result based on the slope bounds of the nonlinearity which was conjectured by Kalman. 


8.6 Examples 


Several examples are presented in this section to allow comparisons of the various approaches. 
Example 1 
Consider a system with an open loop transfer function of G(s) = k(1 — s)/2(1 +s) . Substituting s = jw, with k = 1, gives 


_.  Iil-jo _ Ud-je)’  Ud-ow) _— jo 
CGO) = 201+jo)” 204+0) 204+) (+o) 


so that the Popov locus is 


2 TAYe és a ws 

G(jw) = or Ta 4 7 - which on writing Ggw) = X + jY and eliminating w from the 

expressions for X and Y gives the straight lineY = X — 0,5. The locus starts at (0.5,0) for w=0 and finishes at (-0.5, -1) 
for w=o . Figure 8.5 shows the Nyquist and Popov loci, with the latter drawn for positive frequencies only. It is clear that 
a line can be drawn through the origin with slope approximately 2.0 so that the Popov locus lies to its right. This means 
a k equal to infinity, however the criterion also requires lim _, G(s) >— k' , which gives —-2°"' >— k', thatisk <2. 


The Hurwitz sector is (-2, 2) so that the Hurwitz and Popov criteria are the same for positive gains. 
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Figure 8.5 Nyquist and Popov Plots 
Example 2 


Consider a nonlinear feedback system with G(s) =— 1/(1—s7)(1+s%) and a nonlinear element in the sector 
n(x) € (ki, ke). It is required to find conditions for stability of the closed loop system. The characteristic equation for the 
linear system with a loop gain of kis si'Z, + s(7i — LZ) + k — 1 = 0 so that the Hurwitz sector is (1, ©) provided T >T,. For 
the nonlinear system the Popov criterion cannot be applied directly as G(s) is not stable, that is A. & { Ao} for k=0. However 


if a pole transformation is used with a gain p then G,(s) = 1/(° TE 4+s(f—L)+o-1 and a,(x) €[k —0,k: — eo]. 


Now Gp(s) is stable provided T,>T, and p>1 and since the Popov locus of G (jw) is entirely within the lower half plane 
the Popov sector for n, (x) is (0,0c) which gives stability for n(x)€(1,°°) , provided T,>T,, in agreement with the Hurwitz 
sector. For n(x)€(1, ©) the circle D for the circle criterion is a unit circle with centre at (-0.5, 0). It is easy to show that 


the Nyquist plot of G(s) will lie within this circle for all T,, T, so that stability cannot be proved. 


Example 3 


s+0,8 ; 
sG@ 401s + 03549) with n(x) €{N} and {M} such that n(x) € {0,k} and 


n' (x) € {0,k}. The Nyquist and Popov loci for the transfer function are shown in Fig 8.6 from which it can be seen that 


Consider a system with G(s) = 


the Hurwitz sector is greater than the Popov sector. As w — 0, the real part of G(jw) tends to -7.81 so that the locus will 
lie to the right of a vertical straight line passing through (-7.81,0), giving a gain for stability of k<1/7.81=0.128. A Popov 
line drawn to keep the Popov locus to the right intersects the negative real axis at around -3.40 so to satisfy the Popov 
criterion k<0.29. Similarly a line drawn to keep the Nyquist locus to the right will intersect the negative real axis at -0.48 


giving k<2.08 for stability. Finally for the linear system the maximum gain for stability from the Hurwitz sector is 2.61. 
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8.7. Conclusions 


In this chapter some absolute stability results have been presented primarily those using frequency domain criteria. Their 
major advantages are that they allow comparisons with the approximate DF method and because of their formulation 
allow stability to be proven for classes of nonlinearity, rather than a specific one, so that they are robust to changes in the 
nonlinearity. Their disadvantages are that results of their application can prove conservative for a specific nonlinearity 


and they are restricted in their applicability to a feedback loop with a single nonlinear element. 
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9 Design of Nonlinear Control 
Systems 


9.1. Introduction 


In agreement with the situation for books on linear control systems the early chapters of this book have necessarily been 
devoted to analysis methods rather than design. In practice the design process involves consideration of the approach to 
be used to meet the given specifications, the selection of the components to be used and their modelling and finally the 
use of analytical methods and probably simulation to show that the design meets the specifications. This is normally an 
iterative process requiring a return to earlier stages in the process, making changes and then continuing. The difficulty if 
any components are nonlinear is, as has been seen, that the behaviour is dependent on the magnitude and type of input 
signals to which the system is subjected, no necessary and sufficient conditions for stability assessment are known and 


no general system analysis method is available. 


Thus if the system specification requires the input step response to have certain properties, then it is not sufficient to know 
that this will be the case for a step input of a single magnitude, as with a linear system, but one must determine that the 
requirements will be satisfied for all step inputs within the allowable range. Also if the system is subject to multiple inputs 
it is not sufficient to determine the behaviour for the inputs separately since nonlinear interactions between them may 
cause unexpected behaviour. The behaviour of high order nonlinear systems is dependent on the region of the state space 
in which operation takes place, as has been demonstrated for a second order system by use of the phase plane. Once a 
system has been designed simulation provides an excellent tool for checking the performance within the allowable state 
space region of operation, however, if it is done without theoretical understanding it can prove inefficient and possibly 


erroneous if the simulations fail to cover appropriate scenarios. 


One might argue that the first objective of a design should be to try and add compensation to make the system linear, so 
as to allow easy analysis, but this may not be possible or even desirable. The advantage of linearity is that the response 
is independent of the input amplitude but it may only be possible to have linear behaviour between a specific input 
and output. A very common form of nonlinearity is plant actuator saturation, for example motor torque saturation in 
a speed or position control system. It might be possible to avoid this by using a larger motor but if the system satisfies 
the specifications with a smaller motor then the consequence of changing to a larger one is a more costly design. An 
objective of any engineering design is to meet the specifications at the lowest cost. The specifications are paramount and 
as well as performance and reliability issues they may cover other aspects such as operating costs, weight or volume, etc. 
and environmental concerns. In some cases when the actuator to a plant saturates it can cause problems as the feedback 


loop is then essentially operating open loop. One extreme case is when the plant model is an unstable transfer function. 
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In many cases where nonlinearity occurs it may be possible to do a control design based on an approximate linear or 
linearized system model and then check by simulation that the effects of the nonlinearity still allow the required performance 
to be met. Such a design may involve parameter selection for fixed controllers like PID, or one of several state space design 
methods. The latter may use state variables estimated using observers and although unlike for linear systems there is no 
rigorous justification provided by a separation principle the procedure is effective. Today most controllers are implemented 
using microprocessors so that it is just as easy to implement a nonlinear control law as a linear one. When should this be 


done and why? What are the advantages and disadvantages? 


In the following sections some design methods for nonlinear systems are mentioned and hopefully these will go some 
way to answering these questions. The topics, however, are only covered briefly to give a ‘flavour’ of what is available 
and also because to do otherwise would require the introduction of additional theoretical concepts. Complete texts are 
available on most of the topics and references will be given in the appropriate sections. There are also several books that 
give a good overview of many of the topics considered in this chapter, such as references 9.1 and 9.2. The latter, reference 
9.2, because of the author's industrial background is particularly interesting in offering a much more practical perspective 


than is often the case. 


9.2 Linearization 


Nonlinear systems are often designed using linear methods based on an approximate linear model. Simulations may 
show that satisfactory performance may be obtained from a design using one linear model but in some situations this 
will not be possible due to the large changes in the best linear model with respect to the system operating point. Linear 
controllers may then be designed for the different operating point models, or linearizations. The controller parameters 
are then ‘scheduled’ in some manner based on knowledge of the operating point. Many approaches are possible which 
can be confirmed to work satisfactorily in simulations although there is sometimes a lack of rigorous theory to support 


the techniques. 


The approach is particularly common in the aerospace industry where a common approach is to employ state feedback 
with the gains scheduled according to the specific operating point. The design is typically done off line in that a linear 
mathematical model is obtained for each operating point and then the required controller gains calculated and stored 
ready for implementation around that specific operating point. Various methods are also used to avoid discontinuous 
changes of controller gains with change in the operating point. Two good introductory references for work in this field 
are [9.3, 9.4]. Ifa new mathematical model is identified during operation then this provides a standard form of adaptive 
control. All control designs typical involve considerations of system performance under plant changes and how severe 
these may be expected to be will affect whether a robust design can be accomplished by a single controller or whether 


scheduled controllers or an adaptive control is required. 


Many processes which need to be controlled in the chemical industry often have quite nonlinear dynamics so these are 
sometimes controlled by PID controllers which allow the three PID terms to be scheduled according to the operating 
point. In recent years there has also been a significant amount of work on the topic of multiple models for use in both 


modelling and control [9.5] 
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9.3 Frequency Response Shaping 


In classical controller design methods for linear systems it is acommon procedure to add known forms of compensator, such 
as phase lead or lag, to achieve a desired frequency domain property, such as a minimum allowable gain or phase margin. 
In such designs the requirement is then to obtain suitable controller parameters so the system meets the specifications. 
This is again a satisfactory approach for a nonlinear system and in terms of the DF approach, assuming a single nonlinear 
element, becomes one of shaping the open loop frequency response locus so that it lies at an appropriate distance from 
the C(a) locus. The C(a) locus may be regarded as a "moving (-1, 0)’ point on the Nyquist plot so a gain and phase margin 


may be defined for each value of a. 


Typically a design will involve shaping the frequency response locus to ensure that the value of the gain and/or phase 


margin for all a is never less than a given value. 


Apart from the obvious way of doing this by comparing the location of the C(a) and G(jw) loci on a Nyquist plot, or on 
a Bode plot, particularly for a single valued nonlinearity, another approach is to make use of methods for determining 
limit cycles. Using such an approach the amount of gain or phase shift, the latter by adding a time delay, which may be 


added to the loop before a limit cycle results can be calculated. 
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9.4 Nonlinear Compensation 


Given the block diagram for a feedback system with a nonlinear element, possibilities may exist, dependent upon the 
location of the nonlinearity, for compensation using another nonlinearity. For example, if nonlinearity exists at the input 
to the plant then it may be possible to include a nonlinearity in the controller output to make the loop behaviour ‘more 
linear. Typically this might be used to provide a higher gain for small signals to compensate for a valve dead zone. From 
purely block diagram considerations it is possible to include another path to the system output in parallel with the plant, 
but in practice this is not physically possible. However, the closed loop design problem may be made linear by taking 
such a path and connecting it to the feedback signal as shown in Figure 9.1. If in the figure the nonlinear element N, is 
taken equal to 1-N, then it is easy to show that the closed loop characteristic equation is linear as is the transfer function 


from the input, R, to the control signal, U, which is 


U(s) _ ssoG 
R(s) 1+GG. 


(9.1) 


r(t) 


Figure 9.1 A Nonlinear Compensation Method 


Sometimes nonlinearity occurs in a minor feedback path, for example, when the damping is a nonlinear function of speed, 
unlike viscous friction which is assumed to vary linearly with speed. Obviously, if speed is measured then feeding back 
a nonlinear function of this signal can be used to make the total speed feedback linear. Practically, however, this means 
knowing the form of the nonlinear speed damping and may also require an additional sensor for speed measurement or 


the implementation of a speed estimator. 


A very common use of nonlinear compensation is for the prevention of integral windup, which is provided with most 
commercial PID controllers. As mentioned previously saturation always eventually occurs in a plant, whether it be the 
maximum torque from a motor, the maximum power from a heater, the maximum flow through a valve and at this stage 
further increase in the controller output will not increase the level of the system output, allowing of course for any possible 
delay in the response. If the controller contains an integral term then its output will increase in magnitude and not start to 
decrease until its input changes sign. This can cause large overshoots in the system response and in the case of an unstable 
plant, since the loop is basically open, possibly catastrophic behaviour. Thus, integral windup protection is achieved 
by stopping the integration when the controller output signal reaches a certain level. Today, since PID controllers use 
microprocessors various software algorithms are used to implement the procedure, but in the days of analogue controllers 
a simple approach was to incorporate a nonlinear feedback path around the integrator. Another approach to avoiding 
integral windup [9.6] is simply to switch on the integrator only when near to the steady state, which in some ways seems 
more logical as the integration normally has a destabilising effect on the loop and the only reason it is usually included 


is to eliminate steady state errors. A slight difficulty is in obtaining a suitable algorithm to implement the switching on 
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and off conditions without causing transient problems. 


9.5 Compensation using DF Models 


In the previous section it was assumed that any nonlinearity in the control system could be represented in a simple block 
diagram structure. In practice, such modelling may not be easily achieved or may not be of relevance for the approach 
because of the complexity of the plant model structure with multiple dynamic and nonlinear elements. In this case a 
possible approach is to make use of a frequency dependent DF model, N(a,q), for the plant, which may be obtained from 
either a mathematical model of the plant or possibly frequency response identification testing. As long as the frequency 
response shows a dominant output at the input frequency the model should be satisfactory for use in compensator design 
studies. Two approaches are possible, the first simply being a DF implementation of the structure of Figure 9.1, that is if 
the nonlinear plant has the model, Mi(a,w), then the DF model fed back around the controller N2(a,@), should be such 
that Mi(a,@) + No(a,w) + G(jw) a linear model. To DF accuracy the open loop is then essentially linear, so that linear 


compensation can be used. 


An alternative approach is to place a series compensator before the plant which has a DF model N2(a,w) such that the 
series combination N,(a2,@) and Ni(a@,@) are approximately equal to G(jw) a linear model. The different subscripts on 
the a’s have been used to denote that the amplitudes into the DF models will not be the same. More details on the actual 
approach to achieving this design can be found in reference 9 7. The approach has also been extended for use with relay 
autotuning (see chapter 7) to implement nonlinear terms in a PID controller. When the relay autotuning approach is used 
for a nonlinear plant then the estimated critical point will change as the relay output amplitude is changed. This not only 
enables one to estimate the ‘degree of nonlinearity’ in the plant but also enables the design of a nonlinear PID to produce 


roughly a linear input output response [9.8]. 


9.6 Time Optimal Control 


Early work in time optimal control was done for second order systems by making use of the phase plane approach. 
However, the topic developed rapidly in the 50’s and a large number of papers on the subject were published in the 60’s. 
Although much of the early work was done in the USSR, major developments started in the USA in the 60’s primarily 
driven by the space program where an important requirement was to develop control strategies using minimal fuel, a 
problem which is closely related to time optimal control. The concepts used are based on classical Hamiltonian mechanics 
and generally are known as Pontriagin’s maximum principle or the Bellman-Hamilton-Jacobi approach. The books by 
McCausland [9.9], Jacobs [9.10], and Glad and Ljung [9.1] give good introductions and more advanced books are those 
by Ryan [9.11], Athans and Falb [9.12] and Lewis [9.13]. The mathematical theory of these approaches is not easy so here 
a simple introduction to the topic is provided by way of the phase plane for a double integrator plant. The solution to the 


problem using the classical Hamiltonian approach is given in chapter 6 of reference 9.9. 
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When considering the initial condition response of a second order system in section 2.3.2 it was seen that the phase plane 
plot for a closed loop system consisting of a relay and a double integrator plant transfer function consisted of parabolas. 
This corresponds physically to a position control system, assuming a constant inertial load and no friction, being driven 


at constant acceleration, K, with the equation of the parabolas being 
XxX} — X35) = 2K (x; — X10) (9.2) 


where x, is position and x, velocity. If the input to the relay is — x, — Ax. the switching between parabolas corresponding 


to positive and negative acceleration (deceleration) takes place on the line —x, — Ax, = 0 in the phase plane. 


It is easily appreciated that the fastest possible way to the origin from an initial condition of (-h, 0) is to have the maximum 
acceleration of +K to a distance -h/2 from the origin followed by the maximum deceleration of -K to the origin. The 
final trajectory to the origin is the parabola x; = — 2Km, or if the motion starts from (h, 0), rather than (-h, 0), the final 
trajectory is x7 = 2Kx. Thus to achieve a minimum time response in the ideal case these two parts of parabolas should 
form the switching curve not the line —x, — Ax, = 0. A minimum time response will therefore be obtained if the input 
to the relay is —x2|x.|— 2Kx,. This has been called a V|V| system [9.14]. It is, however, often implemented by having 
the switching boundary as — x. + /2K sgn(—™),/|—m| = 0, when it is known as the SERME (sign error root modulus 


error) system. 


The optimal solution is seen to be dependent on knowing the system model, in this case the parameter K assuming 


negligible friction. If K is estimated, however, the switching boundary can be adjusted by some adaptive strategy. 
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9.7. Sliding Mode Control 


The phenomenon of sliding was encountered in Example 2 of Chapter 2, where a phase plane plot was constructed for a 
second order system containing a relay with dead zone and hysteresis. What happens, as explained in the example, is that 
the dynamics of the motion changes when the relay switches and it is possible that at both sides of the switching line the 
resultant dynamic motion is towards the switching line, producing a motion which is directed with rapid switching along 
the switching line. The interesting feature recognised about this sliding motion is that once the sliding motion is started 
the trajectory is defined and may not change even with some changes in the system dynamics. Sliding mode control is 
a particular form of variable structure control, where the feedback law is selected from several possibilities by decision 
rules. Variable structure control systems (VSCS) were first studied in the USSR [9.15]. A second order system will again 


be considered in the next section before a few comments on the more general situation. 


9.7.1 A Second Order System 


Consider again the same system as Example 2, Chapter 2, but in this case with the relay nonlinearity replaced by an 
ideal relay. Again the feedback to the relay is assumed to be —x — Ax: and the transfer function of the plant after 
the relay to be K/s*. The Simulink diagram is shown in Figure 9.2 with A = K = 1. Assuming an initial condition of 
(— x, 0) then since the relay output wu = sgn(— x — Ax2) the initial value of u is positive and the equation of motion 
is x3 = 2K(% — Xo) which meets the switching line x + Ax = 0 at say (1,2). The motion after switching is then 
described by x3 — 22, =— 2K(x, — x) and this will move back to the switching line if its slope is less than — 1/A. That 
is if — K/x2 < — 1/A, which gives 


x2 < AK. (9.3) 


This condition for ensuring sliding is known as the reachability condition. 


Sane 


Transfer Fen Scope2 


Integrator 


Gain2 Scope3 


Figure 9.2 Simulink diagram for a system with a sliding mode. 


It is normally given in a more general form since if the input to the relay is taken as s so that s = — x, — Ax) = 0 is the 


switching line, the control signal, u, is given by 
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u = sgen(s) (9.4) 


then if the switching line is to be reached it requires lim,.o.s < 0 and lim,.o-s > 0, which can be expressed as the single 


reachability condition 


ss <0 (9.5) 


As the motion in the phase plane is symmetrical the actual reachability condition using equation (9.3), or as can easily 


be shown from equation (9.5), is 


When ideal sliding takes place the resulting motion is along the switching line and is governed by the differential equation 


om + ax = 0. The control action required to bring about this motion is 


discontinuous and not defined on the switching line. However, in addition when sliding takes place s = 0 and therefore 


s = 0, which in the case of the example is 
5s =—% — AKu = 0, giving 

u=—%/KA =— x%/AK (9.7) 
which is known as the equivalent control. 
Figure 9.3 shows the phase plane response for the system of Figure 9.2 from the initial condition (-1, 0). The switching 
line x: + %2 = 0 is first met with x < 1, the reachability condition of equation (9.3) for A = K = 1, so sliding takes place. 
Figure 9.4 shows the system output, x, (Scope ); — x2 (Scope3) and the approximate equivalent control (Scope 2) obtained 


by filtering, u , which is shown in Figure 9.5 continuously switching during sliding. The approximate equivalent control 


can be seen to become equal to —%; = — x. during sliding as given by equation (9.7) for A = K = 1. 


oa 
-1 D5 5 O27 26 25 4 43 2 1 i) 


Figure 9.3 Motion in the phase plane for input from (-1,0) 
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Figure 9.4 Signals in the Simulation 


Figure 9.5 Control signal u 


Let it now be assumed that the plant model is not an ideal double integrator but has some nonlinear damping. This is 


illustrated in Figure 9.6 where the nonlinear damping has been made quite significant, whereas in practice one would not 


expect such a significant change from the assumed model, so that the change in the response on the phase plane is clearly 


seen in Figure 9.7. The switching line is reached on a different trajectory but once it is attained sliding takes place as before. 
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Figure 9.7 Phase plane response for system of Figure 9.6. 
The reachability condition for this example, from equation (9.6), is |x2|< 1 and it may be that reachability is required 
for a larger region in the phase space. Consider the control signal u to be a summation of state feedback and the signum 
function, so that 
Uu = kx + kx + osgn(s) 


Then ss = S(— xX — Ax) = S(— x AK[ kim + ko + sgn(s)]) 


Choosing ki = 0 and ky =— 1/AK gives ss = s[—AKosgn(s)] =—AKo| s | =—7|s| where 7 is a positive design scalar. 


Selection of ss < — 7|s| is known as the 7-reachability condition. The control signal 


U = — (%/AK) + osgn(s) (9.8) 
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which is seen to be the equivalent control of equation (9.7) plus the signum function. In this case the switching surface is 
reachable from all regions in the phase plane. Figure 9.8 shows the Simulink diagram for this system and Figure 9.9 the 
response from an initial condition of (-6,0). Curve (a) is for the control law of equation (9.4), that is the system of Figure 
9.2 and curve (b) for the control law of equation (9.8), that is the system of Figure 9.8. Also shown for the second case 


is curve (c) the response from (-6, 2). 


= 


Scope 


XY Graph 


Figure 9.8 System with control law of equation (9.8) 


4 5 =~ 3 2 -1 o 1 


Figure 9.9 Phase plane responses for the two different control laws. 


From the considerations of a simple second order system it has been seen that when the reachability condition, given by 
equation 9.5, is satisfied, sliding is obtained down the switching curve, which of course need not be linear. Further the 
motion can be regarded as having an equivalent control action to produce a sliding response which is first order, that is 


reduced by one, from that of the original system. 


The formulation presented has used the ideal mathematical signum function, which cannot be realised in practice, although 
it can be nearly ideal in some implementations. However, the effects of very rapid switching may cause practical problems 
such as rapid wear of components, so the ideal may be approximated by using hysteresis in the switching or using an 


approximation to the signum function such as sgn(s) = where 0 is a small positive constant. 


|s|+0 
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9.7.2 Further Comments 


The presentation of the previous section has hopefully indicated the objectives of sliding mode control design, namely to 
design a controller so that the motion will reach, and subsequently remain on a sliding surface, in contrast to a line for 
a second order system, in the state plane. The dynamics of the system when confined to the surface is described as an 
ideal sliding motion and for a system with a single input will be reduced in order by one from that of the original system. 
The controller design for high order systems is not easy to cover in a short presentation of the topic but there are many 


specialised texts available on the subject such as Utkin [9.16],Zinober [9.17], and Edwards and Spurgeon [9.18]. 


It is shown in these texts how designs can be achieved using various state space approaches for MIMO as well as SISO 
systems, what type of plant uncertainty will not affect the sliding motion, and how the ideas can be applied to nonlinear 
systems and sliding surfaces, not just hyperplanes. Like most state space approaches all the states may have to be available 
for implementation, so additional instrumentation or observers may be required, although some results are available 
based only on output feedback. There has also been significant coverage in recent years on the design of observers using 


sliding mode principles[9.18]. 
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9.8 Fuzzy Logic 


The concept of fuzzy logic was introduced by Lotfi Zadeh [9.19]. It is perhaps an unfortunate choice of name from a control 
design viewpoint because a controller designed using fuzzy logic is deterministic having a defined output for a defined 
input. It is the procedure of producing the output for a specific input which is the ‘fuzzy aspect. This procedure varies 
slightly according to the approach used but that due to Mamdani [9.20] contains membership functions, fuzzification 
and defuzzification procedures, and a set of rules. Typically, standard forms for the first three may be available and the 
designer has only to provide a suitable set of rules. Thus, as stated by Zadeh, it enables a person with no knowledge of 
control theory to design a controller by writing a few logic type rules, which in principle may be selected based on the 


observed response of the plant. 


The first point to make about the procedure is that in general the controller designed will be nonlinear, so the wisdom of 
doing this if the plant is essentially linear is debatable. Secondly, it seems to be a fairly common practice to write the logic 
rules for a controller in the error channel in terms of the error and its derivative. Thus the controller becomes a nonlinear 
controller of the form f(e,é). It is not easy to see what advantage a controller of this PD form may have compared with the 
less general PD form fi(e) + £(é), a classical PD controller with nonlinear P and D terms. Many papers describing a fuzzy 


logic controller end up with a nonlinear one yet fail to check the response variation with the amplitude of the input signal. 


A very good discussion on fuzzy logic controllers can be found in the results of the DynLAB project of reference 9.21. In 
particular it explains that the output is simply a static nonlinear function of the input(s) and how membership functions 
affect the shape of the resulting nonlinearity. Finally there seem to be no examples given in the literature of obtaining 
fuzzy logic controllers for plant transfer functions which are known to be difficult to control, for example, plants with 
right half plane zeros, resonant modes or which are unstable. Fuzzy logic ideas have been used in other areas of science 


and engineering; a text covering fuzzy logic in control is reference 9.22. 


9.9 Neural Networks 


Neural network techniques provide a mechanism for obtaining a static or dynamic nonlinear structure using various 
types of nonlinear and linear elements, and configurations. Typically parameters in the structure are changed adaptively 
(training) to try and match given input and corresponding output data structures. Probably the earliest use of this concept 
in control was a variant on this theme where the objective was to produce a nonlinear function, thus adjusting the system 
parameters to produce a known output for a given input. This was done to achieve required nonlinear functions on analogue 


computers where typically the individual functions were a nonlinear element consisting of two adjustable linear segments. 


Using sets of input and the corresponding output data, training techniques are used in a neural network to adjust the 
network coefficients, nonlinearity in the static case, so that they model the system producing the data sets. To test the 
resulting model it is usual to subject it to new sets of input data to check how closely the resulting output matches the 
corresponding known output data from the given system. Obviously the success of the approach, as we have seen from 
the interesting behaviour of nonlinear systems discussed in chapter 1, will be highly dependent on the structure used 
in the network relative to that of the actual system. This will affect many factors such as the accuracy of the model, the 


required training time and indeed whether the behaviour can be modelled, for example in the case of chaotic motion. 
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Although this type of black box approach can be useful, the importance of being able to incorporate physical knowledge 
into system modelling, particularly the structural connections between linear and nonlinear elements, cannot be over 
emphasised. The simple example of illustrating the approximation of a static nonlinearity by a neural network provides 
a good illustration of this. Suppose one has two symmetrical odd nonlinear elements (i) ideal dead zone and (ii) a 
cubic. Further assume there is a choice of neural network structures for doing the modelling consisting of either (a) 
linear segmented nonlinearities or (b) power law characteristics. Both of these functions may be described as infinite 
approximators, as they can approximate any nonlinear function exactly by taking an infinite number of them. Thus if 
the approximators (b) are used for (i) and (a) for (ii), then indeed to get a good approximation a large number will be 
required. However, if (a) is used for (i) and (b) for (ii) then only a small finite number will be required to get an exact 
result. For example, if the power laws are x" with n = 1, 3, 5 etc only the cubic term is required for (ii). Neural networks 
have been used successfully in many areas of science and engineering so there are many books on the topic. Two dealing 


with their application in control engineering are references 9.23 and 9.24. 


9.10 Exact Linearization 


In section 9.3 the possibility of making a system linear by the addition of other nonlinear elements was discussed from 
a block diagram viewpoint. Some of the situations mentioned are special cases of an approach, based on state space 
methods, known as exact linearization. The topic is considered in depth in Isidori [9.25] but introductions are given in 
several textbooks including, Slotine and Li [9.26], Corriou [9.27] and Glad and Ljung [9.1]. This section will follow quite 


closely the presentation given in the last reference. 


The basic concept of exact linearization is to provide a feedback to the plant input, normally taken as u, which makes the 


response from the new input, r, linear. Consider the state space description 


x1 = X2 
Xo = M(x) + (x) + u (9.9) 
ye mK 


If r is a reference input then a closed loop control can be taken with 


“= r—m(m) — mim) — GX — xX (9.10) 


which gives the closed loop system 


Xi =X 
X= AX —-hxy+r (9.11) 
y=ram 


This is a linear system between the input r and output y and the values chosen for a, and a, allow location of the closed 
loop poles. When the state space description (9.9) has y as linear position the technique is often referred to as ‘computed 


force’ or when y is a rotary position as ‘computed torque’ and is often used in robotics. 


On the other hand if the equations are 
X= X2 + mx) 


X. =—X—X.+U (9.12) 


yeHmu 
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Then it is not possible to directly linearize n,(x,) by the choice of u since the two quantities are not in the same equation. 
Although if the block diagram is drawn it is easily seen how the system can be linearized if x,can be measured or estimated 
using an observer and an additional input given to x,. x,, however, is typically an internal variable which cannot be 
influenced directly. It is easy to see that a general form of state equation which will allow exact feedback linearization if 


all the states are available, known as the controllability canonical form, with state vector x, is :- 


Xi X2 
Xp X3 
= : yeu (9.13) 
Xp F(x) + b(x)u 
with the input u chosen to be 
u = (1/b(x))(r — kx — f(x)) (9.14) 


Here u consists of three terms, the new input r, the linear state feedback kx and a term to cancel the nonlinearity. The 


last line in the state equation becomes 
Xn = r— kx (9.15) 
So the equation is linear between r and y and the state feedback gain vector k can be chosen to suitably place the system 


poles. This looks all right in principle but one obvious difficulty is if for some x, b(x) becomes small or zero. To consider the 


situation further a more general form of nonlinear state space representation will be considered and some definitions given. 
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9.10.1 Relative Degree 


Consider a system described by the general differential equation 


x = f(x) + ug(x) 
y = h(x) 


(9.16) 


where f, g and h are infinitely differentiable functions of the nth order state vector x. Obviously y does not depend explicitly 


on u so that any changes in y will occur gradually via x. 
Differentiating the output equation one obtains 
y = h,(x)x = hy (x) (f(x) + g(x) u) (9.17) 


where /,(x) denotes the partial derivative of h(x) with respect to x, that is 


oh oh oh 
h,(x) = ( ae (9.18) 


Equation (9.11) can be written 

y=h.f + uhg (9.19) 
where the x dependence of h., f and g has been omitted for simplicity. It can be seen that y now depends directly on u 
for at least some values of x, if and only if h.g # 0, and in this case the system will be said to have relative degree 1. If 
h,g = 0 then one can differentiate again to obtain 


y = (Af) xf + uf) (9.20) 


The system is said to have relative degree 2 if h.g = 0 and (h,f).g 4 0 so that » depends directly on u. Otherwise one 


can continue differentiating until a derivative of y is directly dependent on u. 
It is useful to introduce some partial differential operators for the repeated differentiations, namely 


—¢_9 Oo = ae 
L; =f ae + jeeethige ax, and Ly = 81 Ae =F ecsiscaves 8n ae (9.21) 


L, is called the Lie derivative in the direction of f. If v denotes the relative degree, which is the lowest derivative of y that 
depends on v, then from the above results it can be seen that 
v=1:L,h 40 
v=2:Lh=0 Lhe£0 (9.22) 
v=3Lbh=0 L2hA40 


etc. 


One thus has the following definition 
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Definition 1 


The relative degree v is the least positive integer such that L,L'h #0. 


It can be shown that this definition is consistent with that for a linear system where the relative degree is the difference 


in the degree between the polynomials of the numerator and denominator terms of the transfer function. 


9.10.2 Input-Output Linearization 


From the above for a system of relative degree v, one has 
y? = Dh +ul,Ue'h (9.23) 
so that if the feedback, assuming L,L; 'h # 0, 
1 y 
is introduced then one has 


yur (9.25) 


a linear relation ship between the input r and the vth derivative of the output. Here Lh and L,L';'h are functions of x 


and in particular if the latter is zero for any value of the state x the required control signal becomes infinite. 
If the coordinate change 
Z = (x) (9.26) 


is made with 


by h 
go Lyh 
& |_| Lyth 


(9.27) 
Dyer arbitrary 


pn arbitrary 
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Then the new system description will be of the form 


v4 22 

Zz z3 

Zz | _ |é() + un(2) (9.28) 
Za Wi (z,u) 

Zi Wn» (Z, U) 

y Z 


for some functions &, 7 and . The linearizing feedback is 


u=(r—€)/n (9.29) 


which gives the closed loop system 


ral Fé) 
va @ 
a ae r (9.30) 


Zo th(z,(r — €)/n) 


Zn Yn» (Z, (r— €)/n) 
y rar 
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Thus the whole system has not been linearised and there is possibly some nonlinear dynamics affecting the states 
Civ igor: ,Z,. This dynamics is not visible in the output and is called the zero dynamics of the system. If the above 
calculations are performed for a linear system then it is possible to show that the zero dynamics are the zeros of the 


transfer function from u to y. 


9.10.3 Exact State Linearization 


It will be noted from the above that the zero dynamics disappear when v = n. In this case there is a state feedback which 


makes the system exactly linear in the state variables 
ZHVi B= Vues eh y' 


and the resultant system will be in the controllability companion form given in equation (9.7) in z, that is 
ca £2 


— : (9.31) 
Zz, E(z) + un (z) 


y Zi 


and can be linearised with the feedback u = (r — &)/7 and the n' line of the above equation is z, = r. 


9.10.4Examples 


Two examples are now considered to illustrate the method. 


Example 1 


Consider the nonlinear system 


x = n(x) + % 


2. =—-M+U 


Obviously the input- output linearisation is dependent on the output y which affects h in equation (9.27). First let the 


output be taken as 
y=Hum+X. 


Differentiation of y shows that the first derivative is dependent on u so that the system is of relative degree 1. From equation 


(9.27) this means y = z = x + x. and z can be chosen arbitrarily. Choosing z, = x gives 
Z=n(m)+u=n(e)+u 


B=n= n(x) +X. = n(z) +4 -2 
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Choosing u = r —n(z), gives Z, = r and the nonlinear zero dynamics is given by the above equation for z, and will be 


stable if z, is bounded. 


Now as a second choice consider the output to be y = x then differentiating gives 
y=x =n(m%) +x and 
Y=n' (mq) me +X = n' (m1) (N(%1) +2) — 2 + U 


where n'(x) = ae Thus choosing u = r+ x. — x2.n'(m) — n'(m)n(x) gives y = r. In terms of z 
a=X 

=X = N(x) +» = Z and 

Z=r 


if the feedback u is taken as 
u=r—an'(ay+a—n(a). 
More directly from equations (9.24) and (9.27) 


Z=y= and 


&=Lh=f hs +f Ohe = (m(x) +) = % 


OX 
Lite ln +x] ah Ono +] =0+1= land 
Lih =f lnc + 2] +f tO a = [n(x) oF x|n' (x) — X2 giving 


U=Pr+X —xn'(%) — n'(x:) n(x) as before. 


Consider now that it is desired to place both poles for the closed loop response to z, at -1 then choosing 
u=r—7y-—2m-—Zzn'(4a)+z—Nn(z) gives 
pat), Dey 

~\-1 -2 1 


which has the required characteristic equation of s* + 25+ 1=0 
In particular if n(x) =— x? then the feedback in terms of z is 
u=r1—%7—244+3n7+a4+U 


and Figure 9.10 shows the simulink diagram for the system. The feedback in the bottom part of the diagram corresponds 
to the last three terms in the above equation for u, which without further terms would make z, = r. There are two positive 
feedback nonlinear terms in these three to make the system linear between the input r and z, and the other two terms are 
the linear ones to place both poles at -1. The linear responses to z (Scope) and z (Scopel) for a step input as expected 


correspond to those of transfer functions 1/(s + 1)* and s/(s + 1)’, respectively and that to z is shown in Figure 9.11. 
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Figure 9.10 Simulink diagram for example 1 
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Figure 9.11 Input-output step response for example 1. 


This appears to be a rather complicated but interesting way of obtaining a linear input-output response but what happens 
if there are modelling errors? To illustrate the problems for this particular case let the actual nonlinearity be n(x) = — ax? 
with 10% variations in a, that is with a equal to first 0.9 and then 1.1. Step responses are shown in Figure 9.12 for step 
magnitudes of 1 and 1.14, respectively. The four responses are marked in the boxes with the first number being, a, and the 
second the step magnitude. Although three responses are satisfactory the fourth for a step input magnitude of 1.14 with 
a = 0.9 actually becomes unstable. This is not surprising in view of the strong positive feedback, whereas the stability of 


the original system is not a problem. 


Figure 9.12 Step responses for Example 1 with nonlinearity parameter changes. 
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Example 2 


Consider the third order system 
X =X. — ax 

Xo = — bx. + n(xs) 

X3 =— cx, + U 


y=m 


where (x3) is any nonlinear element. 
This means 


Xo — aX 0) 
h=m f= [-s + vs g= [| so that 


— CX3 1 


Lh = fi =X. — ax 
Lih = f(—a)t+ ff =-— am + am — bx + n(x) 


Lyh = f(a’) + f(— (a + b)) + fn) 
=—ax+(a@ + ab+ b’)x. — (at b)n xs) — cx3n' (x3) 


and 
L,L-h = n' (xs) 


Thus the required feedback to make z, = r is 


1 


u=— {r+ ax —(@ + ab + b’)x + (a + b)n(xs) + cxsn' (xs) } 
n (x3) 


or in terms of z apart from 71' (x3) 


= 


= n' (x3) 


{r + aba + (a + b)z + cxsn' (x3)} 
Thus if the three poles for the linear system in z are all to be placed at -1 then the required u is given by 
u= ao + abe + (a+ b)% 4+ cx3n' (x3) — Z -— 3m — 3z}. 


Here if n(x3) = x3 then n'(x3) needs to be limited to a small value when x; = 0 to limit the value of the signal u. It is also 
clear that the method can easily be used if the nonlinearity is linear segmented as its slope is known. Figure 9.13 shows 
the implementation in Simulink for the case of a = 0 and b = c = 1 and the linear segmented nonlinearity 


n(xX3) = Xs for | x3|< 1 
= 0.3x; + 0.7 for x3 > 1 
= 0.3x; — 0.7 for x3<-—1 


The Scope blocks are connected to the components of z and the output of the product block is the signal u. The relay 


elements connected in parallel are used to give n' (xs). 
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Figure 9.13 Simulink implementation for Example 2 


9.10.5 Further Comments 


To summarise, the procedure described in the previous sections for a system with relative degree, v, equal to the dimension 
of the state x allows a nonlinear transformation to be found from x to a new state z and a control signal u which will make 
z linear with respect to the new input r. A feedback controller can then be designed in terms of the original state x, or z, 
to control the linear dynamics of z with respect to r. If the relative degree, v, is less than n the dynamics of some states 
of z, known as the zero dynamics will remain nonlinear and need to be checked in case they are unstable. The required 
control signal u is given in equation (9.24) and since its denominator is a function of the state x it may not be finite at 
all points in the state space. The method requires the state x to be available so these values must either be measured or 


estimated using a nonlinear observer. 


Perhaps the major advantage of the method is that it allows a stable design of a nonlinear system but requires an accurate 
mathematical model. When the assumed model does not accurately match the system it may be difficult to assess the 
effects it may have and as shown by the first example the effects may be undesirable. Also as mentioned previously it may 
not be appropriate to remove the effects of the nonlinearity, in this case by using control energy, as their existence may 


not be detrimental to the system performance. 
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9.11 General Comments 


Several methods for the design of controllers for nonlinear systems have been briefly covered in this chapter. There is a 
huge literature on many of the topics so a few appropriate references have been cited but thousands more can be found 
by putting topics and key words mentioned in the above into Google. Comparison of the various approaches is not easy 
except perhaps in the case of optimal control where a solution is obtained which meets a certain criterion, and therefore 
by definition is optimal for that case. However, even in this case, the criterion may have had to be perturbed from the 
true requirements in order to obtain a mathematical solution or, alternatively, it might be known that the model used 
has deficiencies. A large number of published papers are usually devoted to a particular approach not to comparisons of 
methods whilst in practice the approach used will be influenced by many factors not all of which will be easy to describe 


analytically. 


Finally, one must again stress the importance of simulation to support studies of nonlinear control systems no matter 


what analytical theories are used. 
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